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ABSTRACT 


V 


r  r 

The  aspect  of  '•Robustness*  for  linear  multivariable  systems  is 

^  7 

analyzed  in  time  domain.  Both  'Stability  Robustness'  and  'Performance 

f 

Robustness*  are  combinedly  considered  to  meet  stability  and  performance 
requirements.  First  a  stability  robustness  condition  in  time  domain 
(in  terms  of  eigenvalues)  is  presented  and  examples  are  given  which  in¬ 
dicate  that  the  proposed  robustness  condition  is  less  conservative  then 
the  corresponding  frequency  domain  condition  as  well  as  another  recently 
proposed  time  domain  condition,  both  given  in  terms  of  singular  values. 

V 

Next  a  technique  is  presented  to  further  reduce  the  •conservatism'of  the 
proposed  condition.  A  design  algorithm  that  incorporates  both  *stability 

y  y  77 

robustness'  and  'performance  robustness ^  into  the  design  procedure  suggested 
ip  the  summer  faculty  program  report ,  is  modified  with  the  help  of  new 
definitions  of  robustness  indices.  Computer  software  to  implement  the 
algorithm  is  presented  along  with  simple  examples  to  illustrate  the 
concepts.  Based  on  the  experience  gained  by  the  minigrant  research,  areas 
of  future  research  are  recommended. 


Nomenclature 


Real  Vector  Space  of  Dimension  a 
Dirac  Delta 

Belongs  to 

Frequency  Variable 

Spectral  radius  of  the  matrix  [.] 

The  largest  of  the  modulus  of  the  eigenvalues  of 

[.] 

Singular  values  of  the  matrix  [.] 

Eigenvalues  of  the  matrix  [.  ] 

Symmetric  part  of  a  matrix  [.] 

Modulus  Matrix  *  Matrix  with  modulus  entries 
for  all  i 


I.  INTRODUCTION  AND  OBJECTIVES 


It  is  well  known  that  the  inaccuracies  in  the  mathematical  models  of 
physical  systems  can  severely  compromise  the  resulting  control  designs. 

The  errors  associated  with  mathematical  models  of  physical  systems  may  be 
broadly  categorized  as  i)  parameter  errors,  ii)  truncated  models  (errors 
in  model  order),  iii)  neglected  or  incorrectly  modeled  external  distur¬ 
bances  and  iv)  neglected  nonlinearities.  It  is  the  inevitable  presence 
of  these  errors  in  the  model  used  for  design  that  eventually  limits  the 
performance  attainable  from  the  control  system  designs  produced  by  either 
classical  (frequency  domain)  or  modem  (time  domain)  control  theory.  The 
problem  of  model  errors  is  more  critical,  in  general,  for  large  scale 
Linear  Quadratic  Gaussian  (LQG)  optimal  control  problems  and  in  particular, 
for  aerospace  applications  like  Large  Space  Structure  (LSS)  control  [1] 
and  other  aeroelastic  systems  [2],  These  applications  are,  of  course,  of 
extreme  importance  to  the  U.S.  Air  Force.  The  fundamental  problem  of  these 
Distributed  Parameter  Systems  (DPS)  control  is  the  control  of  a  large  di¬ 
mensional  system  with  a  controller  of  much  smaller  dimension  (model /control¬ 
ler  truncation)  compounded  with  modal  data  uncertainty  (parameter  errors) . 

In  the  light  of  these  observations,  it  is  evident  that  ’robustness'  is  an 
extremely  desirable  (sometimes,  necessary)  feature  of  any  feedback  control 
design  proposed  for  DPS  control.  'Robustness'  studies  of  Large  Scale  LQG 
regulators  is  the  central  theme,  of  the  present  research. 

For  our  present  purposes  a  'robust'  control  design  is  that  design 
which  behaves  in  an  'acceptable'  fashion  (i.e.  satisfactorily  meets  the 
system  specifications)  even  in  the  presence  of  modeling  errors.  Since  the 


system  specifications  could  be  in  terms  of  stability  and/or  performance 
(regulation,  time  response,  etc.)  we  can  conceive  two  types  of  robustness, 
namely,  'Stability  Robustness'  and  'Performance  Robustness'.  Limiting  our 
attention  in  this  research  to  'parameter  errors'  and  'model/controller 
truncation'  as  the  two  types  of  modeling  errors  that  may  cause  instability 
(or  performance  degradation)  in  the  system,  we  fornally  define  'stability 
robustness'  and  performance  robustness'  as  follows: 

'Stability  Robustness':  Maintaining  closed  loop  system  stability  in  the 
presence  of  modeling  errors  mainly  parameter  variations  and  model/controller 
truncation. 

'Performance  Robustness':  Maintaining  satisfactory  level  of  performance  in 
the  presence  of  modeling  errors  mainly  parameter  variations  and  model/control¬ 
ler  truncation. 

Implicit  in  the  definition  of  'Performance  Robustness'  is  the  require¬ 
ment  of  'stability'  for  performance  robustness  studies.  However,  performance 
robustness  studies  which  require  (or  assume)  stability  may  only  have  limited 
application.  On  the  other  hand  concentrating  design  efforts  on  'stability' 
alone  is  not  prudent  because  the  performance  requirements  may  not  be  met  by 
that  design.  Thus,  simultaneous  consideration  of  stability  and  performance 
in  the  design  process  is  more  appropriate.  Most  of  the  current  published 
literature  addresses  either  the  '  stability  robustness'  aspect  or  the  'per¬ 
formance  robustness'  robustness'  aspect  separately.  Most  of  the  interesting 
work  on  'stability  robustness'  is  done  in  frequency  domain  using  singular 
value  decomposition  [3-11]  while  much  of  the  useful  research  on  'performance 

robustness'  is  carried  out  in  time  domain  using  sensitivity  approaches  [12-14]. 


Design  studies  that  treated  the  stability  robustness  aspect  in  time  domain 
and  studies  which  combined  both  stability  robustness  and  performance  ro¬ 
bustness  into  the  design  process  have  been  scarce.  Towards  this  direction, 
research  was  initiated  on  these  aspects  by  the  author  during  the  Summer 
Faculty  Research  Program  period  (Summer  ’82)  and  consequently,  a  stability 
robustness  condition  in  time  domain  and  a  design  algorithim  that  incorporates 
both  stability  robustness  and  performace  robustness  into  the  design  process 
were  proposed  [IS].  Later  as  part  of  Mini  grant  work  the  following  research 
was  proposed. 

i)  To  probe  further  into  the  possible  refinement  of  the  stability 
condition. 

ii)  To  develop  computer  software  for  automating  the  design  algorithim 
and  finally 

iii)  To.  illustrate  the  methodology  by  examples  representative  of  Large 
Space  Structure  models. 

In  what  follows,  a  summary  of  the  research  that  accomplished  the  above 
tasks  is  presented  followed  by  a  discussion  of  areas  of  further  research. 

II.  Work  Done  During  the  Mini  Grant  Period 

IIA.  Time  Domain  Analysis  of  Stability  Robustness 

Much  of  the  published  literature  treats  stability  robustness  of  feed¬ 
back  control  systems  in  the  frequency  domain  with  the  help  of  singular 
value  decomposition  [3-11  §  16].  While  the  singular  value  analysis  is  a 
useful  tool  to  generalize  the  Nyquist  criterion  for  the  multi-variable  case, 
it  is  not  the  only  way  to  characterize  the  stability  of  a  perturbed  system, 


particularly  for  systems  described  by  state  space  models.  For  systems 
described  by  state  space  differential  equations,  the  analysis  of  stability 
in  the  time  domain  becomes  a  viable  alternative.  The  time  domain  ro¬ 
bustness  analysis  may  have  a  number  of  advantages  over  frequency  domain 
treatment  as  explained  in  later  sections.  The  time  domain  treatment  is 
more  or  less  analogous  to  the  frequency  domain  treatment  in  spirit,  but 
with  the  stability  conditions  given  in  terms  of  eigenvalues  rather  than 
singular  values  since  eigenvalue  analysis  is  more  appropriate  for  time 
domain  stability  assessment.  As  stated  earlier,  it  is  of  interest  to  note 
that  eigenvalue  analysis  was  used  even  in  the  frequency  domain  treatment 
(Ref.  [16]). 

In  what  follows,  we  first  present  the  main  mathematical  result  (as  a 
new  theorem)  which  forms  the  basis  for  developing  a  condition  for  the 
stability  of  a  perturbed  matrix.  This  result  is  then  extended  to  the  case 
of  a  linear  system.  These  results  are  further  specialized  to  the  case  of 
LQG  regulators  with  perturbations  in  the  form  of  i)  parameter  variations 
and  ii)  truncated  models. 

Main  Result: 

Let  F  and  E  be  two  real  square  matrices. 

Theorem  1:  If  F  is  negative  definite,  then  the  matrix  F  +  E  is 
negative  definite  if 

-1 

P  {  IW  ]s  }  <  1  (1) 

Proof:  Given  in  Appendix  A. 

The  above  theorem  has  an  interesting  implication.  It  suggests  that  if 


a  given  matrix  is  written  as  the  sum  of  a  negative  definite  matrix  F 
and  a  perturbation  matrix  E 


i.e.,  F^  =  F  +  E  where  F  is  negative  definite 
then  by  use  of  theorem  1  one  can  get  a  condition  for  the  negative  de¬ 
finiteness  and  hence  the  stability  of  matrix  F^ . 


slication  to  Linear  State  Space  Models: 


Let  us  consider  the  linear  time  invariant  system 

n 

X  =  A  x  ,  x  (0)  =  xQ  ,  x  £  R  (2) 

where  the  'nominal*  matrix  A  is  asymptotically  stable.  Let  there  be  an 
additive  perturbation  E  in  the  matrix  A  so  that  the  perturbed  system 
matrix  is  A  +  E.  We  are  now  interested  in  the  stability  of  the  matrix 
A+E.  The  straightforward  application  of  theorem  1  gives  the  following 
result . 

Theorem  A:  The  perturbed  system  matrix  A+E  is  stable  if 


-1 

p  [(EnAd  >.J  <  1  (3) 

where  A  =  A.  ♦  A  A,  *  Diag  [a.],  i=l,2...n,a.  is  any  real  negative 

d  G  Q  X  X 

entry  (a.<0)  and  E  =  (A  +  E)  . 

n  v  e  's 

Proof:  Since  A  is  asymptotically  stable,  we  can  always  write 


A  *  A ,  +  A 
d  e 

where  A^  =  Diag  [c^J ,  i  =  l,2,...n.,  is  real  and<0. 

thus  A  +  E  *  A,  +  (A  +  E) 
a  e 


-  5  - 


(4) 


where  is  a  symmetric  negative  definite  matrix.  Applying  theorem  1 
gives  the  result  of  theorem  A  since  a  negative  definite  matrix  is  always 
a  stability  matrix. 

At  this  stage,  it  is  appropriate  to  comment  on  the  applicability  of 
these  stability  conditions.  It  is  to  be  noted  that  the  derived  stability 
condition  is  not  particularly  useful  when  one  knows  both  the  matrices  A 
and  the  perturbation  matrix  E,  in  which  case  the  stability  of  the  matrix 
A  +  E  =  A^  +  Ag  +  E  is  determined  by  simply  looking  at  its  eigenvalues. 
However,  in  a  practical  situation,  one  doesn't  exactly  know  E.  One  may  only 
have  knowledge  of  the  magnitude  of  the  maximum  deviation  that  can  be  ex¬ 
pected  in  the  entries  of  A.  In  that  case  the  entries  of  E  are  such  that 


(CD 

where  is  the  magnitude  of  the  maximum  deviation. 

The  following  theorem  enhances  the  usefulness  of  the  proposed  stability 


criterion. 


Theorem  B:  The  perturbed  system  matrix  A+E  is  stable  for  all  perturbations 
E^  satisfying  (Cl)  if 

P  [(EmAm  )s]<  1  CS) 

where  A  =  A^  +  A^,  A^  =  Diag  [ou],  i  =  1 ,2, . .  .n.o^c  0 


E.  -  dAJ  *  ws  \  ■  lAJ 

Proof:  Observe  that 

-1  -1  -1 
PCf  (|Ae|  ♦  A)s  |Ad|  }s]  >  p[|{  (Ae+  E)s(Ad)  >s|]  >  p[{  (Ae+  E)^)  }g] 

(for  all  E„  satisfying  (Cl)).  (7) 


6 


This  follows  from  the  fact  that 


i)  for  any  given  square  matrix  F 


P  C|F|)  >  p  (F)  (Ref.  [16]) 


(8) 


ii)  for  two  given  square  non-negative  matrices  F  and  F^  such 
that  Fuj  >  F2..  for  all  i.j 

p  (Fx)  >  p  (F2)  (Ref.  [17]) 


(9) 


thus 


0  (AOs'  <  1  *  »  t(E„Ad‘l  ),J  <  1 


(10) 


which  in  turn  implies  A+E  is  stable  for  all  satisfying  (Cl) . 

Q.E.D. 


Thus,  in  practice,  we  test  the  condition  of  theorem  B  and  if  satisfied 
we  can  guarantee  asymptotic  stability  for  all  perturbations  E_  satisfying 
(Cl). 


Note  that  there  is  some  flexibility  in  splitting  the  stable  matrix  A 

into  A,  and  A  .  One  immediate  choice  for  A ,  could  be 
a  e  d 

Ad  *  Re  [Ai (A) ] ,  i  =  1,2, - n  (C2) 

We  now  compare  the  ability  of  the  proposed  time  domain  condition  in 
predicting  the  stability  of  a  perturbed  system  with  i)  the  corresponding 
frequency  domain  condition  and  ii)  another  time  domain  condition  recently 
proposed,  both  given  in  terms  of  singular  values. 

Recently,  in  Ref.  [5],  Lee  et  al .  considered  the  state  feedback 
regulator  of  Fig.  1  and  provided  a  condition  for  the  stability  of  an 
additively  perturbed  system.  It  is  shown  that  for  a  system  which  is 


initially  stable,  stability  will  be  ..rintained  as  perturbations  A  A  are 
added,  provided  that 


Thus  the  frequency  domain  condition  fails  to  predict  the  stability  of  the 


perturbed  system  matrix  A  +  AA. 

Li  Li 

Proposed  Time  Domain  Condition: 


A 


d 


and  thus  E  = 
m 


so 


0  1.5 

and  A  = 

8  0 

m 

1.5  0 

0  1 

J 

-1 


[(E  A 
1  m  m 


>,J  ■ 


0  0.84375 


0.84375  0 


=  0.84375  <  1 


(15) 


Thus  the  perturbed  system  matrix  A  +  AA  is  guaranteed  to  be  stable  and  it 

LL 

is  indeed  stable. 

Comparison  with  the  time  domain  condition  of  Lee  [18]: 


In  [18],  Lee  proposes  the  following  stability  robustness  condition  in 
time  domain  in  terms  of  the  singular  values .  The  perturbed  system  matrix 
A  ♦  AA  is  stable  if 

a  (AA)  <  -a  .  (A)  cos  (©  .  ) 
max  '  min'  '  '  mm' 

where  0  .  is  the  smallest  principal  phase  of  A  measured  counterclockwise 
mxn  c 

from  the  positive  real  axis.  The  above  theorem  is  stated  under  the  assumption 
that  the  orthogonal  matrix  U  in  the  polar  decomposition  of  A,  viz 


A  =  U  Hr  or  A  =  Hl  U 


is  a  stable  matrix.  Specializing  this  condition  for  a  symmetric  stable 
matrix,  we  get  the  condition  as 


Example  2:  Let  us  consider  the  same  example  as  before.  The  time  domain 
condition  of  (12)  yields 


a  [AA]  =  1.5  *  I  A. (A) |  .  =1 

max1  J  r  1  i  J  'rain 

Thus  the  above  condition  fails  to  predict  the  stability  whereas  (as  shown  in 
Example  1)  the  proposed  time  domain  condition  succeeds. 

Also  note  that  the  proposed  time  domain  condition  of  this  research 
doesn't  require  the  assumption  of  the  orthogonal  matrix  U  being  stable  which 
involves  testing  of  yet  another  condition  [18]. 


Improvement  of  'Optimism*  of  the  proposed  condition: 


The  flexibility  in  splitting  the  stable  matrix  A  into  A,  and  A  can  be 

a  e 

utilized  to  further  improve  the  optimism  of  the  proposed  condition. 

One  suggested  technique  is  as  follows: 

Procedure :  1)  Write  A,  =  -a I 


where  a>o  is  a  scalar  variable. 
The  stability  criterion  then  becomes 


where  E®  *  d\J+  I E  i  v 

The  above  condition  is  then  tested  for  various  values  of  a  which 
clearly  increases  the  probability  of  guaranteeing  the  stability  of  a  per¬ 
turbed  system  whenever  it  is  stable. 


-8  0 

0  2 

Example  3:  Let  A^  = 

0  -1 

;  AA  = 

2  0 

*  m 

J 

10 


with  A,  =  Diag  [A. (A-.)],  the  proposed  time  domain  condition  becomes 

Q  1 


0.25  0 


1.125 


1.125  0 


=  1.125  {  1 


But  with  A^  =  -al,  the  test  matrix  for  the  proposed  time  domain  condition 


becomes 


P  a 


|a-8 1  2/a  *j 


choosing  a  =  8,  we  get 


0  0.25 


p  *  0.9414  <  1  -*•  the  system  is  stable. 

0.25  0.875 

Thus,  the  proposed  time  domain  condition  with  'a  splitting'  is  able  to 
predict  the  stability  of  the  perturbed  system  and  thus  yields  a  less  con¬ 
servative  result. 

Evidently,  there  is  much  scope  to  improve  the  'optimism'  of  the  proposed 
time  domain  condition  by  an  appropriate  selection  of  the  A^  matrix  and  more 
research  is  warranted  in  this  direction. 


Extension  to  LOG  Regulators 


We  now  extend  the  above  analysis  to  the  case  of  large  scale  LQG  re¬ 
gulators  having  i)  parameter  variations  and  ii)  truncated  modes  as  the 
modeling  errors  (or  perturbations) .  The  strategy  adopted  is  to  model  these 
perturbations  as  additive  perturbations  to  a  nominally  stable  matrix  and  then 
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apply  the  above  eigenvalue  analysis  to  arrive  at  the  stability  conditions. 


Let  us  consider  a  continuous  linear  time  invariant  system  described  by 


*  (t)  =  A  x  (t)  +  B  u  (t)  +  D  w  (t)  ,  x(0)  =  x 


(16a) 


y  (t)  =  C  x  (t) 


(16b) 


z  (t)  =  M  x  (t)  +  v  (t) 


(16c) 


where  the  state  vector  x  is  nxl,  the  control  u  is  mxl,  the  external 


disturbance  w  is  qxl,  the  output  y  (the  variables  we  wish  to  control) 


is  kxl  and  the  measurement  vector  z  is  £xl.  Accordingly,  the  matrix  A 


is  of  dimension  nxn,  B  is  nxm,  D  is  nxq,  C  is  kxn  and  M  is  JUn.  The 


initial  condition  x(0)  is  assumed  to  be  a  zero-mean,  gaussian  random 


vector  with  variance  X  ,  i.e. 

o 


E[x(0) ]  =  0,  E[x(0)  x  1  (0)]  =  X 


Similarly,  the  process  noise  w  (t)  and  the  measurement  noise  v  (t) 


are  assumed  to  be  zero-mean  white  noise  processes  with  gaussian 


distributions  having  constant  covariances  W  and  V  respectively,  i.e. 


E[w(t)]  =  E[v(t)]  *  0 


-  i:  - 


TV,1. 


■  k’-  -  • 


Let  the  above  system  be  evaluated  for  any  control  u  by  the  quadratic 
performance  index 

J  =  lim  -j-  E  [(yTCO  Q  y  (T)  mT(x)PcR0  u  (t)  ]  dx  (20) 


where  scalar  p_  >  0  and  Q,  R  are  (kxk)  and  (mxm)  symmetric,  positive 
c  o 

definitive  matrices,  respectively. 


For  the  case  of  a  deterministic  system,  the  following  modifications 
in  the  system  description  are  in  order: 
i)  Dw  *  0,  v  *  0 

T 

ii)  the  initial  condition,  x(0)  =  X0»X0X0  =  x0 

and  the  index  J  of  (20)  reads 

J  *  (J/'a’[yT(t)  Q  y(t)  +  uT  (t)  pcRq  u(t)]dt  (21) 


If  the  state  x(t)  of  the  stochastic  system  is  estimated  as  a 
function  of  the  measurements  we  assume  the  state  estimator  to  be 
of  the  following  structure 


I 


x(t)  =  A  x(t)  +  G  z(t) 


(22) 


I 

r 

r 

r 

i 

■’ll 

5: : 

L 

[ 
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where 

z(t)  =  z(t)  -  M  x(t) 

is  called  the  'measurement  residual'.  For  a  'minimum  variance' 
requirement,  the  estimator  of  (22)  is  the  standard  Kalman  filter 
[19].  We  refer  to  the  system  presented  in  this  section  as  the  'Basic 
System' . 

We  now  use  theorem  B  as  the  basis  for  extending  the  stability 
conditions  to  the  case  of  LQG  regulators  for  various  cases  of  modeling 
errors  involving  parameter  variations  and  model/controller  trunca¬ 
tion.  These  are  discussed  in  [15].  For  brevity  we  consider  two  cases 
here. 

Case  1:  Parameter  variations  alone,  no  model/controller  truncation: 

For  this  case,  the  following  assumptions  are  made  with  respect 
to  the  model  described  by  equations  (16). 

Assumption  1:  The  matrix  pairs  [A  ,  B]  and  [A  ,  D]  are  completely 
controllable  and  the  pairs  [A  ,  C]  and  [A  ,  M]  are  completely 
observable. 


L 

L 


Since  there  is  no  model/controller  truncation  the  full  order 
optimal  control  for  nominal  values  of  the  parameters  is  given  by 


u  =  G  x  =  - 


A 

K  x 


(23a) 
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V 


m 


X=Ax+Bu+G(z-Mx),  x(0)  =  0 


«  (A  ♦  B  G  •  G  M  )  x  +  G  : 


1  T  -1 
G  =  —PM  V 

Pe  0 


and  P  and  K  satisfy  the  algebraic  matrix  Riccati  equations 


KA  +  ATK  -  KB  —  BTK  ♦  CTQC  =  0 
pc 


T  IT»  * 

PA  ♦  AP  -  PM  — —  MP  *  DWD  =  0 
Pe 


The  nominal  closed  loop  system  is  given  by 


where  A  =  A  +  BG  -  G  M  and  the  closed- loop  system  is  asymptotically 


stable. 


We  are  now  interested  in  examining  the  stability  robustness  of  the 


closed-loop  system  in  the  presence  of  parameter  variations  alone. 


Let  AA,  AB ,  AC,  AM  and  AD  be  the  maximum  modulus  perturbations  in  the 
system  matrices.  A,  B,  C,  M  and  D  respectively.  Then  the  perturbed 


'V -V 


system  can  be  written  as  (Note:  the  filter  parameters  A^  and  G  and  th 
control  gain  G  are  not  subjected  to  variations  because  they  are 
specified  by  the  designer  as  a  function  of  the  nominal  values  of  the 
parameters)  . 


Ax 

A* 


BG 


GM  A  .0  0 

c  I 


AA  ABG  i+AA  (B+AB)Gj 
GAM  0  p (M+AM)  Ac 


Ax 

A 

Ax 


D  0 

0  G 


AD  0 
0  0 


(25) 


By  application  of  theorem  B  of  section  II,  we  obtain  the 
following  design  observation. 

Design  Observation  1:  The  perturbed  LQG  regulator  system  is  stable  for 
all  perturbations  in  A,B,C,M  S  D  in  the  sense  of  (Cl)  if 


-1, 


p[tE»  \  V  •=  1 


(26) 


where 


Ad  =  Diag  [Real  ^(A^jj  ,  ^  = 


BG 


GM  A 


ACL  =  Ad  +  Ae  ’  Em  =  ' Ae '  +E’  E  = 


CJ  , 
[AA  AB  j  G  | 


g|am  o 


!and  v|Ad!  (27) 


Note  that  in  this  case  both  the  nominally  stable  matrix  and  the 


perturbation  matrix  E  are  functions  of  controller  gains  G  and  G. 


Case_2:  Model  reduction  alone,  no  controller  reduction  and 


no  parameter  variations 

For  this  case,  we  treat  the  model  given  by  (16)  as  the  evaluation 
model.  We  assume  the  order  of  the  model,  n,  to  be  too  high  for  the  control 
u  to  be  determined  and  that  there  is  a  control  design  model  of  dimension 
nR  <  n  given  by 

*r  =  Vr  *  V  *  5r  "•  V*  R 


yR  =  CR*R  <28> 

»r  ■  Vr  *  v 


where  the  above  control  design  model  is  obtained  either  by  a  direct 
truncation  of  the  full  order  model  given  by 


*rt 

XR 

br' 

r  -i 

dr 

u  + 

XT 

.  J 

bt 

dt 

(29a) 


y  =  [CR  CjJ 


z  * 


[Mr  Hj.] 


’XR 


n. 


;  VR 


(29b) 


(29c) 
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or  by  a  partial  realization  of  x  involving  some  model  reduction 
technique  (e.g,  [14]). 


Let  the  full  order  control  for  the  reduced  order  model  be  given  by 
u  =  Gr  5Fr  [30a 


*R  =  AR*R  +  BRU  +  GR(Z  " 


(30b 


*  *R  rR  *  GR  Z  “here  AR  =  *R  *  BRGR  *  Vr 

such  that  the  closed-loop  system  matrix  for  the  control  design  model 
given  by 


ar  brgr 


(31) 


cL 


is  asymptotically  stable.  These  controller  gains  GRand  GR  could  be 
optimal  or  non-optimal  with  respect  to  the  model  (28). 

The  closed-loop  system  for  the  evaluation  model  is  obtained  by 
forcing  the  evaluation  model  with  the  controller  of  the  control  design 
model.  Thus,  we  have 


J 


The  stability  of  the  above  closed-loop  system  matrix  which  we  denote 


as  is  to  be  established. 

At  this  juncture  we  assume  the  matrix  Ap  to  be  an  asymptotically 
stable  matrix  (which  is  a  reasonable  assumption  for  large  space  structure 
models).  In  order  to  derive  the  condition  for  stability  of  the  closed- 


loop  system  matrix  ACL?, 


we  write  ACL2 


as 


where  G  and  G  are  the  control  gain  and  estimator  gain  matrices,  re- 
R  K 

spectively .obtained  by  using  the  reduced  order  model  (A  B  C  M  D  }, 

K  j  K  y  K  I  K  y  K  y 

and  are  such  that  the  resulting  closed  loop  system  matrix  (of  the  re¬ 
duced  order  system  which,  of  course,  is  the  first  partition  of  Ag  )  is 

A 

asymptotically  stable.  One  choice  of  G  and  G  could  be  the  standard 

K  K 

optimal  control  gains  of  the  reduced  order  model  {  A  B  C  1  D  } 

K  y  K  y  K  |  K  y  K  y 

under  appropriate  conditions.  Note  that  Ag  is  thus  a  stable  matrix 
and  that  £p  L  basically  contains  the  terms  that  cause  spillover. 


Design  Observation  2:  The  closed-loop  system  matrix  A  is  stable  for 
all  control  design  models  GR  and  GR  such  that 


GRij  I  * 1  W|  GMj| 

if 

13  "Em  Os1  <  1 

where  Ad  +  Ag  -  Ag  c>  Ad  =  Diag  [Re  (As  C)J 
Eb  .  (|Ael  *  6E)S  .  Am  =  |Adl 


and  AE  = 


0  |2BrGr| 

laVV  |2*R! 
IS-rI  IbtgrI 


I A 


RT1 


(34) 


(35) 
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Measure  of  Stability  Robustness 


It  may  be  noted  that  a  variety  of  controller  gains  may  satisfy  the 
proposed  stability  condition  for  given  perturbations.  In  order  to  compare 
different  models/controllers  from  stability  robustness  point  of  view,  it  is 
clear  that  there  is  a  need  for  some  form  of  a  measure  of  stability  robust¬ 
ness.  To  this  end,  we  now  define  a  measure  of  stability  robustness  called 

"Stability  Robustness  Index",  Bc  D 

O  •  K* 

6S.R.  4  I  I  Re  Xmin^ACL^'l  R®  Xmin(ACLP^  ^  j  lReXmin^ACL^  ( 

where  is  the  nominal  closed  loop  system  matrix  and  A^p  is  the  perturbed 
closed  loop  system  matrix  (i.e.  the  closed  loop  system  matrix  formed  by  the 
perturbed  matrices  of  plant,  estimator,  etc.);  here  it  is  assumed  that  the 
controller  gains  are  such  that  the  condition  for  stability  is  satisfied  and 
thus  the  perturbed  closed  loop  system  matrix  is  stable.  The  motivation  be¬ 
hind  this  definition  is  that  the  index  is  a  measure  of  the  magnitude  of  the 
deviation  in  the  minimum  eigenvalue  modulus  of  the  perturbed  system  from  the 
nominal  system. 

By  this  difinition  8C  D=  0  corresponds  to  a  highly  robust  system  from 
stability  point  of  view.  However,  one  aspect  of  further  research  in  this 
development  is  to  investigate  what  is  the  worst  case  deviation  (i.e.  Max 
Bc  D)  for  a  given  set  of  maximum  modulus  perturbations  expected  in  the 
system  matrices.  Once  the  worst  case  ^  is  determined  for  each  controller/ 
model,  say  3-  D  u,  then  this  index  can  be  used  for  comparison  purposes. 

O  •  K  •  W  i 
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Measure  of  Performance  Robustness 


In  similar  lines,  we  define  a  measure  of  performance  robustness  called 


"Performance  Robustness  Index  SD  D  "  i.e. 


BP.R.  -  I  Jp  ■  JN  I / JN 


where  Jp  is  the  value  of  the  performance  index  for  the  perturbed  systen 
(with  A^p  as  the  plant  matrix).  Thus,  Bp  ^  =0  corresponds  to  a  highly 


robust  system  from  performance  point  of  view.  As  mentioned  above,  an  area 


of  future  research  is  to  investigate  the  worst  case  Q  D  for  each  controller/ 

r  •  R  • 


model  and  then  use  this  index  Bp  p  ^  as  the  basis  for  comparison. 


Robust  Control  Desis 


Once  measures  of  stability  robustness  and  performance  robustness  are 


developed,  the  idea  of  a  robust  control  design  is  to  pick  a  controller  that 


gives  a  reasonable  or  satisfactory  trade  off  between  stability  and  performance. 


Specifically,  the  design  algorithms  involves  determining  the  indices  Bc 


and  Bp  «  for  different  values  of  the  design  parameters  p  and  p  which  in 
*  •  c  e 


turn  determine  the  control  and  estimator  gains.  This  information  about  the 

ic 

indices  Bg  p  and  Bp  p  along  with  the  corresponding  nominal  costs  (  T  ) 


and  ^T  can  be  used  to  pick  a  specific  controller  (control  and  estimator 


gain  combination)  as  the  one  which  provides  a  satisfactory  trade  off  between 


stability  and  performance.  The  algorithm  is  thus  iterative  in  nature.  The 


computation  basically  involves  the  use  of  matrix  Riccati  and  Liapunor  so- 


lutions  and  the  eigenvalue  analysis  for  which  standard  easy-to-use  computer 


programs  are  available.  The  details  of  the  algorithms  (in  principle)  are 


given  the  SFRP  report  and  for  reasons  of  brevity,  are  not  reported  here. 


However,  a  brief  account  of  the  "flow-chart"  is  included  in  the  section 


dealing  with  the  computer  software. 
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Discussion  of  the  Theoretical  Development  of  this  Research: 


Some  discussion  about  the  implications  of  the  current  theoretical  develop¬ 
ment  of  this  research  is  now  in  order.  First,  it  may  be  noted  that  the 
proposed  stability  conditions  are  similar,  conceptually,  to  the  frequency 
domain  results  reported  in  Ref.  [1].  However,  there  are  also  some  interest¬ 
ing  differences  between  these  two  (frequency  domain  and  time  domain)  versions. 
Some  preliminary  observations  are  presented  in  the  following  sections. 
Secondly,  these  design  observations  are  useful  in  many  ways  in  both  the  analys 
and  synthesis  of  robust  controllers.  These  are  discussed  in  later  sections. 

Comparison  and  Contrast  Between  Frequency  Domain  Analysis  and  the 
Time  Domain  Analysis 

The  main  differences  between  the  frequency  domain  treatment  and  the 
time  domain  treatment  are  as  follows: 

i)  The  main  differences  between  the  frequency  domain  treatment  and  the 
the  calculation  of  singular  values  of  a  complex  matrix  at  various 
frequencies.  In  the  stability  conditions  of  time  domain,  no  time 
dependence  is  present.  Only  the  eigenvalues  of  a  real  symmetric  matrix 
are  to  be  computed. 


ii)  In  the  case  of  frequency  domain  results,  the  perturbations  are 
mainly  viewed  in  terms  of  'gain'  and  'phase'  changes  [6,7].  In  the 
proposed  time  domain  analysis  the  perturbations  are  viewed  as  'system 
parameter  variations'  and  'system  model/controller  order'  with  con¬ 
stant,  fixed  gains.  It  may  be  noted  that  in  the  time  domain  treatment 
the  nominally  stable  closed-loop  matrix  and  the  perturbed  closed-loop 
matrix  are  both  functions  of  the  constant  controller  gains. 


iii)  In  the  frequency  domain  treatment,  the  requirement  of 
square  matrices  necessitates  the  assumption  that  the  number  of  control 
inputs  be  equal  to  the  number  of  outputs,  which  can  be  satisfied  by 
appropriate  selection  of  the  break  point  of  the  loop.  However,  in  the 
present  time  domain  analysis  no  such  assumption  is  needed. 

iv)  Most  of  the  work  on  robustness  in  the  frequency  domain  con¬ 
centrated  on  the  analysis  problem.  That  is,  analyzing  a  given  closed- 
loop  control  system  to  determine  the  uncertainty  bounds  that  will  make 
the  closed-loop  system  unstable.  This  can  be  done  very  easily  in  the 
frequency  domain.  The  difficulty,  however,  is  that  the  uncertainty 
bounds  can  be  obtained  only  in  the  frequency  domain  and  are  very  dif¬ 
ficult  to  translate  back  into  the  time  domain  to  determine  the  allowable 
perturbations  of  physical  parameters  of  the  system.  This  difficulty 
can  be  eliminated  by  posing  equations  (26)  and  (27)  in  the  framework  of 
a  robustness  analysis  problem. 

Very  little  work  has  been  done  on  the  robustness  synthesis  problem 
in  the  frequency  domain.  That  is,  to  design  a  closed-loop  contorl  system 
to  meet  the  performance  specifications  in  the  face  of  a  specified  set  of 
uncertainty  bounds.  The  proposed  approach  in  the  time  domain  handles  the 
robustness  synthesis  problem  very  well. 

V)  In  the  proposed  stability  conditions,  provision  is  made  for 
considering  different  reduced  order  models  for  control  design  purposes 
(reflected  by  the  presence  of  A^,  B^  matrices  which  could  be  different 
from  A^,B^. . .matrices) .  Just  as  in  the  frequency  domain  different  control 
design  methods  could  be  compared  using  the  singular  value  plots  of  the 
return-difference  matrices,  this  provision  helps  to  compare  different  re¬ 
duced  order  models  (used  for  control  design)  from  a  stability  robustness 
point  of  view. 


vi)  In  the  frequency  domain  treatment,  considering  an  uncertainty, 
for  example,  as  an  additive  perturbation  several  stability  robustness  con¬ 
ditions  can  be  written  which  do  not  imply  each  other  for  practical  systems 
[21].  In  the  present  time  domain  approach  such  difficulty  is  not  present. 
The  perturbations  can  only  be  modelled  as  additive  perturbations  and  yield 
only  one  robustness  test. 

vii)  All  the  norm  bounded  robustness  criteria  in  the  frequency  domain 
are  inverently  conservative  because  they  assume  that  the  phases  of  all  the 
elements  in  the  perturbation  matrix  are  in  worst  possible  direction  which  is 
a  mathematical  extreme.  Some  possible  alternatives  to  reduce  this  con¬ 
servatism  in  the  frequency  domain  are  developed  in  references  [22]  and  [23] . 
The  proposed  time  domain  robustness  criteria  are  also  conservative.  The 
conservatism  enters  the  development  of  the  main  result  because  of  the  re¬ 
quirement  of  the  negative  definitess  of  the  perturbed  system.  The  degree  of 
conservatism,  however,  is  a  subject  for  future  research. 

These  are  some  of  the  preliminary  observations  made  with  respect  to 
the  frequency  domain  and  time  domain  approaches  for  'stability  robustness'. 
Evidently  further  in-roads  have  to  be  made  in  the  investigation  of  this 
relationship  and  this  is  suggested  as  a  future -research  topic.  In  the 
following  section,  the  usefulness  of  the  proposed  design  observations  is 
briefly  discussed. 

Usefulness  of  the  Design  Observations 


The  proposed  design  observations  are  helpful  in  many  ways.  First,  if 
a  nominally  stable  closed-loop  system  matrix  is  specified  (i.e.  either 
or  A  depending  on  the  specific  case  as  discussed  in  section  III),  then 
one  can  use  these  tests  to  determine  the  tolerable  perturbations  in  the 


system  matrices  A,  B  and  M  and  the  order  of  the  modei/controller  before  the 
closed  loop  system  becomes  unstable. 

Conversely,  given  the  perturbations  AA,  AB  ind  AM  and  the  order  of  the 
modei/controller,  one  may  determine  the  controller  gains  to  achieve 
stability  robustness. 

As  indicated  in  the  previous  section,  these  tests  can  be  used  to  conpare 
different  model  reduction  and  control  design  schemes  from  a  stability 
robustness  point  of  view. 

Finally,  these  tests  can  find  applications  in  spillover  reduction 
problems  and  sensor/ actuator  location  problems. 


I  IB,  Computer  Software 

Introduction  r 

The  computer  program  has  been  developed  in  two  packages.  The  first 
package  tests  the  stability  condition  developed  in  the  report  for  two  types 
of  perturbations  viz.  Parameter  Variation  and  Model  truncation.  The 
second  package  then  computes  the  Regulation,  control  and  total  costs  for 
the  open  loop.  Nominal  closed  loop  and  Perturbed  Closed  Loop  Systems. 

The  program  then  computes  the  robustness  indices  SD  D  and  8C  . 

The  program  uses  extensively  the  subroutines  available  in  the  LSLIB 
(Library  for  Control  and  Estimation  of  Linear  Uncertain  Systems)  and  made 
compatible  for  use  in  DEC-10  computer.  The  main  programs  of  the  Stevens 
Computer  Center  and  other  subroutines  are  written  in  Fartran  IV.  The 
LSLIB  package  was  originally  developed  at  Purdue  Aero  Department. 

The  following  paragraphs  outline  briefly  the  program. 

Package  1.  Stability  Evaluation; 

The  algorithm  first  forms  the  nominal  closed  loop  matrix  A^  and  the 
Perturbation  Matrix  AACL* 

a)  Program  PARVA  1  handles  parameter  variation  problem.  The  input 
matrices  are  A,B,C,D,M,Q,R,V,W  and  the  perturbations  DA,DB,DC,DD§DM. 

1)  The  subroutine  SSLQG  solves  the  algebraic  matrix  Rocatti 
equation  for  the  controller  and  estimator  and  gives  the  control  gain  G 

A 

and  estimator  gain  G. 

2)  The  subroutine  MATFTZ  forms  the  nominal  closed  loop  matrix 
and  the  Perturbation  Matrix  AA_. . 
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b)  Program  Trunk  handles  the  Truncated  Mode  problem.  The  input 
matrices  are  AR,  ARJ,  ATR,  BR,  BT>  CRf  CT>  DR,  DT,  MR,  M^,,  Q,  R,  V,  W. 

1)  The  subroutine  SSLQG  solves  the  Ricatti  equation  and  gives  the 

A 

control  gain  GD  and  estimator  gain  G_. 

R  R 

2)  The  subroutine  MATFIT  forms  the  nominal  closed  loop  Matrix  A 

L>  Li 

and  the  Perturbation  matrix  AACL  as  discussed  in  the  report. 

c)  Program  TEST  with  A^  and  AA^  matrices  as  input,  tests  the 
stability  condition  of  the  system.  It  is  interactive  and  permits  one  to 
test  the  stability  of  the  system  for  different  values  of  a  [of  the  a  method 
discussed  in  the  report) .  Once  the  condition  is  satisfied,  the  program 
then  calculates  the  stability-robustness  index  &  of  (  36  )  for  all  the 

o  •  K  • 

controllers  satisfying  the  condition. 

Package  2.  Performance  Evaluation: 

The  input  to  this  program  are  the  matrices  AB,  C,  D,  M,  W,  Q,  DA,  DB, 
DC,  DD,  DM  and  the  scalar  constants  ROEC  and  ROEE. 

1)  The  subroutine  LYAP2  computer  the  open  loop  cost  JQp. 

2)  Subroutine  SSLQG  solves  the  Ricatti  equations  for  the  controller  and 

A 

estimator  and  gives  the  control  gain  G  and  the  estimator  gain  G.  It  also 
computes  the  control,  regulation  and  total  costs  of  the  nominal  closed 
loop  system  (JUN,  JXN,  JN  respectively) . 

3)  Subroutines  MATFT2  and  MATFT4  form  the  matrices  AF11,  AF21,  AF22,  CF, 
DF1 ,  DF2,  and  WF  (Ref.  [15] . 

4)  Subroutines  LYAP  2  and  LYAP  5  solve  the  three  reduced-order  Lyapunov 
equations  (each  of  order  (n  +  n  ))  and  give  the  matrices  PF11,  PF12,  and 


5)  Subroutine  MATFT2  forms  the  matrix  PF. 

6)  Subroutine  MATFT4  forms  the  matrix  QF. 

7)  Subprogram  JRAC2  computes  the  cost  with  the  matrices  PF,  CF  and 
QF  as  input. 

8)  The  algorithm  forms  different  structures  of  the  QF  matrix  each 
one  yielding  a  particular  cost,  thus,  the  control  costs  JUNN,  JUP, 

Regulation  Costs  JXNN,  JXP  and  the  total  costs  JNN  and  JP  for  the  nominal 
closed  loop  and  Perturbed  Closed  Loop  respectively  are  computed. 

9)  The  algorithm  then  computes  the  performance  robustness  Index  BETAPR 
from  the  perturbed  regulation  cost  JXP  and  the  nominal  regulation  cost 
JXNN. 

List  of  Subroutines  Used: 

ABS  -  Makes  all  elements  of  a  matrix  absolute. 

CINV  -  Forms  inverse  of  a  complex  (nonsingular)  matrix. 

DIAMAT  -  Forms  diagonal  matrix  with  real  part  of  eigenvalues  of 

a  given  matrix. 

EIGRF  -  Computes  eigenvalues. 

IDENT  -  Forms  an  identity  matrix. 

T 

LYAP  2  -  Solves  Lyapunov  equation  [form:  AA  *  xx  +  xx  *  AA  +  cc  =  0] 

T 

LYAP  5  -  Solves  Lyapunov  Equation  of  the  [form:  AA  *  xx  +  xx  *  BB 

+  cc  =  0] . 

MADDSB  -  Adds  and  subtracts  Matrices  [form:  A  +  B  -  c] . 

MARB  -  Find  the  biggest  magnitude  of  real  part  of  the 

eigenvalue . 

Find  the  smallest  magnitude  of  real  part  of  the  eigenvalue. 


MARS 


MAT FT  2 

- 

Forms  a  single  matrix  from  4  matrices. 

MAT FT  4 

- 

Forms  a  single  matrix  from  16 

matrices . 

MEQ 

- 

Stores  the  matrix  A  in  the  matrix  B. 

MP  31 

- 

Computes  matrix  product  [form: 

T 

p  =  A  * 

B  *  C 

MP  32 

- 

Computes  matrix  product  [form: 

p  =  A  * 

B  *  CT 

MULRRT 

- 

Computes  matrix  product  [form: 

p  =  A  * 

bt 

MULT 

- 

Computes  matrix  product  [form: 

p  =  A  * 

B] 

SCAMUL 

- 

Multiple  a  matrix  by  a  scalar. 

SSLQG 

- 

Solves  Ricattis  equation. 

STABR 

- 

Forms  and  solves  the  stability 

criterion. 

TRAC  2 

- 

Multiple  2  matrices  and  finds 

the  trace 

of  the 

matrix  product. 

USWCM 

- 

Prints  complex  matrix. 

USWFM 


Prints  real  matrix. 


C  PROGRAM  PARVAt 

C  THIS  PROGRAM  FORMS  CLOSED  LOOP  MATRIX  AND  PETURBATION  MATRIX 
C  FOR  PARAMATER  VARIATION  PROBLEM 

REAL  A<14,14>.8(14.14).C<14.14>.D<14,14>.M<14»14>.DA<14,14) 
REAL  D8<14.14>.DC<14.14).AD<28.28).WU(14.14) 

REAL  DD<14»14>»DM(14»14)*0(14»14)»R<14«14) 

REAL  M<14.14>  .V<  14.14)  tS<  14.14)  .C0C(14.  14  ) 

REAL  K(14»14)»P(14«14)«F(14»14)»G<14»14)» UK ( 300 >  » ZERO ( 14*14) 
REAL  FM< 14. 14) .BG< 14»14) *AC< 14. 14 > . AFM<14. 14) . APG( 14 . 14) 

REAL  ACL (28. 28) .AE(28.28> .DBG( 14. 14) 

REAL  FDM(14«14)»ECL<28»28) .EN<28.28) 

REAL  CC(28.28).DDD(28.28).ENS(28.28) 

C 

COMPLEX  EIGA( 14).EIG(14).EE(28.28).EI(28.28> 

COMPLEX  EIGACL(28).EIGN(28) 

C 

IS-14 

IB-28 

NX-2 

NM-1 

NO-1 

NK-2 

NL-2 

C 

C  GIVE  VALUES  OF  A  AND  DA 

2  READ  (15.*H(A<I.J).J>1.NX>.I-1.NX) 

3  READ  ( IS.*) (  <DA< I.J).J-1.NX).I-1.NX) 

C 

C  GIVE  VALUES  OF  BtDB 

S  READ  <15.t)<(B<l.J).J>l.NM).I*l.NX) 

A  READ  <1S.*X(DB(I.J).J-1.NM).I»1.NX> 

C 

C  GIVE  VALUES  OF  CtDC 

8  READ ( IS. t) <<C(I.J).J-1.NX).I-1.NK> 

9  READ( IS. t) ((DC(I.J).J-lfNX).I-l.NK) 

C 

C  GIVE  VALUES  OF  DtDD 

11  READ  (1S.*X(D(I.J).J*1 .NO) . 1-1 .NX) 

121  READ  (15.*)( (DD( I.J).J-1.NQ).I-1.NX) 


GIVE  VALUES  OF  M1DM 

READ  (IS.*) <  <M< 1 . J> . J-l .NX) .1-1 »NL> 
READ  (1S.*X(DM(I.J).J«1.NX>.I-1.NL) 

GIVE  VALUES  OF  O(NK.NK) 

0  IS  SYM.+VE  DEFN  MATRIX 

READ (IS.*) ((Q(I*J)*J-1»NK)»I-1»NK) 

GIVE  VALUES  OF  ROEC 
READ(15.*)R0EC 
DO  101  1*1 .NM 
DO  101  J*1.NM 
R(I.J)-0 
R(I.I)-1 

GIVE  VALUES  OF  ROEE 
READdS.OROEE 
DO  102  I-l.NL 
DO  102  J-l.NL 
V(I.J)-0 
V(I.I)-1 

READ( IS.*) ((V(I.J). J-l. NL). I-l.NL) 


C  CALLS  IMSL  ROUTINE  FOR  SETTING  OUTPUT  IDENTIFIER 
CALL  UGETI0(3,0.35) 

CALL  USUFM(2HAA. 2. A. IS.NX.NX.4) 

CALL  USUFM(2HDA. 2 .DA. IS.NX.NX.4 ) 

CALL  USUFM(2HBB.2.B.IS.NX.NM.4> 

CALL  USUFM(2HDB.2.DB. IS.NX . NM. 4 ) 

CALL  USUFM(2HCC.2.C. IS .NK .NX . 4 ) 

CALL  USMFM ( 2HDC .2.DC.IS.NK.NX.4) 

CALL  USWFM(1HD.1.D.IS.NX,NQ,4) 

CALL  USWFM(2HDD. 2, OD.IS.NX.NO. 4) 

CALL  USMFM ( 1HM. 1.M.IS.NL.NX.4) 

CALL  USUFM(2HDM.2*DM. IS. NL . NX .4 > 

CALL  USMFM ( 1H0. 1.Q.1S.NK.NK.4) 

CALL  USMFM( 1HM. 1 .M. IS. NQ.NQ .4 ) 

C  rORMS  R-ROEC*R 

CALL  SC AMUL ( ROE C .IS.R.IS.R.NM.NM) 

C  FORMS  V-ROEE *V 

CALL  SC AMUL (ROEE.IS.V.IS.V.NL.NL) 

C  FORMS  MM-D*USDT 

CALL  MP32( IS.NX.D. IS.NO.U. IS.NO.NX.D.IS.MU) 
CALL  HP31(IS.NX.C.IS.NK.0> IS.NK .NX.C. IS. COO 


CALL  USUFM ( 1HR. 1 .R.IS.NM.NM.4) 

CALL  USUFN(1HV.1.V»IS.NL.NL.4> 

CALL  EICRF (A.NX.IS.2.EIGA.EE.IB.UK. IER) 

CALL  USUCH(4HEIGA»4»EIGA* IS. NX. 1.4) 

CALL  SSLGG(XS.A.X6.B. IS.UU, IS.M. IS. V. IS.CGC. IS.R .NX.NM.NL 
1  >0.. FALSE.. UK. IS. 8. IS. K. IS.Q.IS.P. XS.F.EIG.EE.RJ.RJX.RJU) 
CALL  USUFM(1HG»1.G. IS.NM.NX.4) 

CALL  USUFM(1HF»1.F.IS.NX.NM,4> 

CALL  USUFMOHK.1.K.IS.NX.NX.4) 

CALL  USUFM ( 1HP* 1.P.IS.NX.NX.4) 

CALL  MULT(B.G.BG»NX.NM»NX.IS.IS.IS) 

CALL  HULT (F.M.FM.NX.NL.NX.  IS » IS.  IS > 

CALL  MATADD<A»BO»ABQ.NX»NX.IS) 

CALL  MATSUB  < A.FM. AFM.NX dS.IS.IS> 

CALL  EIGRF  < ABG.NX. IS. 2 .EIGN.EE . IB . UK . IER > 

CALL  USUCM  ( 14HEXGENS  OF  A+ BG. 14.EIGN. IB. NX . 1 . 4 > 

CALL  EIGRF (AFM.NX. IS. 2.E10N. EE. IB. UK. IER) 

CALL  USUCM  < 14HEIGENS  OF  A-FM. 14 .EIGN. IB. NX . 1 .4 ) 

CALL  MAD0SB(A.8G.FH.NX> AC. IS. IS. IS. IS)  • 

CALL  MATFT2  < A . BG . FM .  AC . NX . NX . NX . NX , ACL . 2*NX 
1.2ANX.IB.XS.XS) 

CALL  USUFM <  3H ACL  »  3  »  ACL . I B . 2  *  NX » 2*N  X . 4 ) 

NACL-24NX 

NECL-NACL 

CALL  ABS(G.NM.NX.IS) 

CALL  ABS(F.NX.NL.IS) 

CALL  MULT (DB.G.DBG.NX.NM.NX.  IS.  IS .  IS > 

CALL  MULTtF. DM. FIIM.NX.NL. NX.  IS.  IS.  IS) 

CALL  MATFT2 ( DA >  DBG . FDM . ZERO >  NX . NX . NX  > NX • ECL . NECL 
XfNECL.IB.IS.XS> 

CALL  USUFM (lHE.Xt ECL . IB . NECL . NECL . 4 > 

WRITE <36. SI ) < (ECL( X . J) . J-l.NECL) . 1-1 .NECL) 

CALL  MAT  ADD  <  ECL  .  ACL  »  ACL .  NECL .  NECL .  I B  > 

CALL  USUFM  <  5H ACL+E . S . ACL .IB. NECL . NECL . 4 ) 

CALL  EIGRF (ACL .NECL .IB >2. EIGACL.EE. IB. UK. IER) 

CALL  USUCM< 14HEXGEN  OF  ACL+E. 14. EIGACL. IB. NECL. 1 .4 ) 


>* 

s* 


PROGRAM  TRUNK 

THIS  PROGRAM  FORMS  CLOSED  LOOP  MATRIX  AND  PETURBATIDN  MATRIX 
FOR  A  TRUNCATED  MODES  PROBLEM 

REAL  AR<14. 14) »BR< X4» 14) »CR< 14. 14) »DR<14» 14>  »MR< 14. 14). AT < 14.14) 
REAL  BTC14. 14) »CT< 14.14) »AD<28.2B).UH<14.14).A<2B.28) 

REAL  DTC14.14) .MT< 14. 14) »0<14»14).R<14»14) »COC( 14. 14) 

REAL  W(14»14)rV<14»14> »S(I4»14) .ART < 14.14) »ATR( 14. 14) 

REAL  K(14»14).P<14»14)»FR<14»I4).GR(14.14) .UKtSOO) »ZERO( 14 .14) 
REAL  PRMR< 14.14) »BRGR< 14.14 ) »AC< 14.14) .RAFM< 14. 14) »RABG< 14. 14) 
REAL  ACL(28.28) .AE(28.28) »BTGR< 14.14) 

REAL  FRMT<14.14) .ECL (28. 28) >EN(2B.2B) 

REAL  CCC28.2B) .DDD(29.28) .ENS (28.2B) 

COMPLEX  EXGA<28) .EIG( 14) .EE <28. 20) .El (28.28) 

•COMPLEX  EXGACL(28).EIGN<28> 

XS-I4 
IB-28 
•  NXR-2 
NXT-1 

NX-NXR+NXT 

NM-1 

NO-1 

NK-3 

NL-1 

BIVE  VALUES  OF  AR.ART.ATR  AND  AT 
READ  ( 13.XX (AR(I.J).J-l .NXR) .1-1 .NXR) 

READ  <13»*X(AT< X . J) . J-l .NXT) . X-l .NXT) 

RCAD< IS.*) < (ART  < I* J) . J-l >NXT ) . 1-1 .NXR) 

RE AD ( IS.*) ( < ATR< I . J> . J-l.NXR) .1*1 .NXT) 


GIVE  VALUES  OF  BRIBT 
READ  <1S»*X<BR<I»J).J-I.NM).I-1.NXR> 
READ  (1S.*X(BT(X.J).J-1.NM).I-1.NXT) 

GIVE  VALUES  OP  CRtCT 
RCAD( 1S.*)((CR(X.J) .J-l.NXR) » 1-1 »NK  > 
READ(  1S.*X(CT<I.J).J-1  .NXT  ).I-l»NK) 


c 

C  GIVE  VALUES  or  MRIMT 

14  READ  < 15»* ) < (HR< I . J) » J-l .NXR) , 1-1 »NL ) 

15  READ  (1S,*H(MT(I, J) . J«1 .NXT > . 1-1 »NL > 
C 

C  GIVE  VALUES  OF  Q(NK.NK) 

C  0  IS  SYM.+VE  DEFN  MATRIX 

16  READ (13.*) ( <Q< I. J) . J-l.NK) .I-l.NK) 

C 

C  GIVE  VALUES  OF  ROEC 

READ(1S.*)R0EC 
DO  101  I-l.NM 
DO  101  J-  1 .  NM 
R< I . J)*0 
101  R(I.I)“1 


C  GIVE  VALUES  OF  ROEE 
READ( 13. * )ROEE 
DO  102  I-l.NL 
DO  102  J-l.NL 
V<I.J)-0 

102  V(I.I)«1 

C 

C  GIVE  VALUES  OF  U(NO.Na> 

READ(1S.*)< (U(I. J) . J»1,NG)»I-1,NQ) 

c 

C  CALLS  IMSL  ROUTINE  FOR  SETTING  OUTPUT  IDENTIFIER 
CALL  UGETI0<3.0.3S> 

CALL  USUFM (2HAR.2. AR. IS.NXR. NXR ► 4 ) 

CALL  USUFM ( 2HAT .2. AT . IS . NXT . NXT . 4 > 

CALL  USUFM  (3HART.3.  ART  .IS.  NXR  .  NXT .  4  > 

CALL  USUFM < 3HATR . 3 . ATR . IS ► NXT . NXR  > 4 > 

CALL  USUFM (2HBR.2.BR. IS. NXR.NM.4) 

CALL  USUFM(2HBT. 2. BT» IS. NXR.NM.4) 

CALL  USUFM(2HCR.2.CR. IS. NK, NXR. 4) 

CALL  USUFM (2HCT . 2.CT . IS.NK.NXT .4) 

CALL  USUFM(2HDR.2. DR. IS. NXR .NO. 4) 

CALL  USUFM (2HDT .2. DT. IS. NXT .NQ.4) 

CALL  USUFM ( 2HMR • 2 . MR . I S . NL . NXR . 4  > 

CALL  USUFM (2HMT .2. MT . IS.NL.NXT.4) 

C  FORMS  R-ROECtR 

CALL  8CAMUL < ROEC .IS.R.IS.R.NM.NM) 

C  FORMS  V-ROEEtV 

CALL  SCAMUL ( ROEE . I S . V . I S . V . NL . NL  > 

CALL  USUFM ( 1HQ. 1 .Q.IS.NK.NK.4) 

CALL  USUFM ( 1HU.1.U.IS.NG.NQ.4) 

C  FORMS  UU»0*U*DT 

CALL  MP32< IS.NXR .DR . IS. NQ.U. IS.NQ.NXR. DR. IS.UU ) 

CALL  MP31 < IS.NXR. CR.IS.NK. 0. IS.NK.NXR.CR. IS.CQC) 

CALL  USUFM (1HR.1.R.IS.NM. NM . 4 ) 

CALL  USUFMdHV.l.V. IS.NL.NL. 4) 

CALL  MATFT2 <  AR . ART . ATR » AT . NXR . NXT . NXR . NXT . 

1  A. NX .NX. IB. IS. IS) 

19  CALL  EIGRF(A.NX> IB.2.EIGA.EE.IB.UK. IER) 

CALL  USUCM(4HEIGA. 4.EIGA. IB.NX. 1 .4) 

20  CALL  SSLOG( IS. AR.IS.BR. IS.UU. IS. MR. IS. V, IS. COC. IS. R. NXR.NM.NL 
1. 0.. FALSE.. UK. IS. S. IS. K. IS. GF. IS. P. IS. FR.EIG.EE. J. JX. JU) 

CALL  USUFM(2HGR.2.GR. IS. NM. NXR. 4) 

CALL  USUFM ( 2HFR . 2 . FR . I S . NXR . NM . 4 ) 

CALL  USUFM (1HK. 1 .K. IS.NXR .NXR. 4 ) 

CALL  USUFMdHP.l.P. IS.NXR. NXR.4) 

22  CALL  MULT<BR.GR.BRGR.NXR.NM.NXR» IS. IS. IS) 

23  CALL  MULTtFR.MR.FRMR.NXR .NL.NXR. IS. IS. IS > 

CALL  MULT(BT.GR.BTGR»NXT .NM. NXR. IS. IS. IS) 

CALL  MULT<FR, MT.FRMT.NXR.NL. NXT. IS. IS. IS) 

CALL  MATADD < AR.BRGR .RABG.NXR .NXR . IS ) 

CALL  MATSUB ( AR . FRMR . RAFM . NXR  r IS . IS . IS ) 

CALL  EIGRF  (RABG.NXR. IS. 2. EIGN. EE . IB. UK. IER ) 

CALL  USUCM  (17HEIGENS  OF  AR+BRGR. 17 .EIGN, IB. NXR . 1 . 4 ) 
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CALL  EIGRF(RAFM.NXR» IS . 2 . E IGN. EE » IB . UK . IER) 

C 

40  CALL  USUCH  < 1 7HEIGENS  OF  AR-FRMR. 17 .E1GN. IB r NXR . 1 . 4 > 

24  CALL  MADDSB( AR . BRGR.FRMR.NXR.AC.IS.IS.IS.IS) 

C 

NACL-2*NXR+NXT 

NECL-NACL 

C  CALL  MATFT3 

25  CALL  MATFIT(AR.BRGR. ZERO. FRMR .AC . ZERO. ZERO. ZERO. AT, 
1  NXR.NXR.NXT.NXR.NXfi.NXT.NACL.NACL.ACL.IS.IS. IS.IB) 

CALL  USWFM( 3H ACL .3. ACL  *  IB  »2*NXR  +  NXT  »2*NXR4NXT » 4 ) 
WRITE(36»51)( ( ACL( I » J>  » J«1 »NACL) » I“1 .NACL ) 

C  FORM  E  MATRIX 
SL-2 


29 

CALL 

SCAMUL ( SL . IS . BRGR . I S  » BRGR . NXR . NXR ) 

30 

CALL 

BCAMUL (SL » IS. FRMR. IS. FRMR. NXR. NXR) 

31 

CALL 

SCAMUL (SL. IS. AC. IS. AC. NXR. NXR) 

32 

CALL 

ABS( BRGR. NXR. NXR. IS) 

33 

CALL 

ABS (FRMR. NXR. NXR. IS) 

34 

CALL 

ABS( AC.NXR .NXR. IS ) 

37 

CALL 

ABS (ART .NXR.NXT .IS) 

38 

CALL 

ABS( ATR.NXT .NXR . IS) 

41 

CALL 

ABS(FRMT. NXR.NXT. IS) 

42 

CALL 

ABS(BTGR.NXT.NXR.IS) 

43 

CALL 

MATFIT( ZERO. BRGR. ART.FRMR.AC.FRMT.ATR.BTGR. ZERO. 

1  NXR .NXR.NXT  .NXR.NXR .NXT.NECL.NECL  »£CL» IS. IS. IS. IB) 

C 


STOP 

END 


C  PROGRAM  COST: 

C  THIS  PROGRAM  COMPUTES  THE  REGULATION  COST  .CONTROL  COST  AND  TOTAL 
C  COST  FOR  THE  OPEN  LOOP. NOMINAL  CLOSED  LOOP  AND  PERTURBED  CLOSED 

C  LOOP  CASES  AND  THE  PERFORMANCE  INDEX. 

c  dimensions: 

REAL  A(9»9) » DA( 9 .9 ) . B <9. 9 ) »0B( 9 .9) .U(9.9)»M(9,9>.DM<9»9>»V(9.9> 
REAL  a(9.9).R(9,9>.SS(9.9>»KK(9,9).G(9,9>.P(9,9).F(9»9>,A12(9.9> 
REAL  A21 (9.9) . A22(9. 9  > » A32<9 .9) . A33<9. 9) . A34  <9.9  >  »A4 1(9.9) 

REAL  A43(9»9) .AF 11(10.18) »AF12(I8.1B>  »  AF21 ( IB. 18) .AF22a8,18) 

REAL  D(9»9)»DD(9»9)»DF1( 18. 18) .DF2( 18. 16)  »C 1 9» 9) .DC(9,9) ,C23<9»9) 
REAL  CF(36.36) .QF(36.36) .UF( 18* 18) »PR1  (18.18) . PR21 (18.18) 

REAL  PR22< 18.18) »PR2< 18.18) .PR31 ( 18.18) »PR32( 18.18) *PR33( 18.18) 
REAL  PR3( 18. 18) . PF 11(18. 18). PFI2(18.I8) .PF22( 18. 18) .PF(36r36> 

REAL  TPF12( 18. 18) > UK (500 ) . ZER0I9.9 ) *UU( 9. 9 ) rCOC ( 36. 36) . DUM < 18.18) 
REAL  POP  <  9 . 9 ) 

REAL  JOP . JXN . JUN . JN . JXP . JUP . JP . JXNN . JUNN . JNN 
COMPLEX  EIG1 (36) >EIG2<36) .MODI (36.36) 

COMPLEX  I MODI (36.36). M0D2 (36.36). IM0D2 (36.36). UKC ( 1 500 > 

IR-9 

IR2-18 

2  IR4-36 
NN-1 
NM-1 
NO-1 
NK-1 
NL-1 

NN2*2»NN 

NN4-44NN 

NOL-NQ+NL 

3  NKM2*2* ( NK+NM  > 

C  INPUT  values: 

8  READ(50.«)((A(I.J),J=1.NN),I-1.NN) 

RE AD (50.*) ( ( DA < I , J) , J«1 , NN  > . I«1 >NN) 

71  READ(S0.*>( (B(I.J).J>1.NM)»I>1.NN) 

READ (50.*) (<DB(I.J).J*1.NM).I"1.NN> 

72  READ(30.*X(C(I.J>.J»1,NN>.I-1.NK> 

READ<50.*)( (DC(I.J).J-l.NN) ,1-l.NK) 

73  READ(50.*)( < D( I . J) , J»1 . NO ) . I*1.NN> 

READ(50.*)((DD<I.J).J»1.N0> .I-l.NN) 

74  READ(50.*)((M(I.J).J-1.NN).I-1.NL) 
READ(S0.*)((DM(I.J).J-1.NN).I-1.NL> 

76  READ<50.*)((U(I.J).J*i,N0) .I-l.NQ) 

77  READ (SO.*) ( (V(I.J). J*1,NL) .I-l.NL) 

READ (50.*)  ROEC 

READ(S0.«)  ROEE 

READ (50.*) ( (0(I.J).J»1,NK> .I-l.NK) 

READ (30.*)  ITTY 

81  CALL  IDENT(R.NM.IR) 

82  CALL  I DENT (V.NL. IR> 

83  CALL  SCAMUL(ROEC. IR.R. IR.R.NM.NM) 

84  CALL  SCAMUL(ROEE.IR.V.IR.V.NL.NL) 

C  OPEN  LOOP  COST 


c 

501  CALL  MP32< IRfNNfDf IR f NO ,U f IR f NO f NN , D , IR r UU ) 

502  CALL  MEQ( AfDUM  fNN  fNNf  IR  f IR2  > 

503  CALL  EIGRF<DUM.NN.IR2»1»£IC2.H0D1fIR4.UK.IER> 

504  CALL  CINV<M0IUfIM0IHfUKfIR4fNN) 

505  CALL  LYAF2(EIG2fIR4fMGD1f I MODI , IR ,UU , IR  fPOP  fNN , WKC  > 

506  CALL  MP31 <IRfNNfCfIRfNKfQfIRfNKfNNfCfIR2  fDUM ) 

J0P=TRAC2 < IR . POP , I R2 . DUM f  NN  f  NN  > 


CALL  MP31(IRfNNfCf IR,NK ,0 f IR fNN fNNfC f IR4 fCQC > 

5  CALL  SSLOQ<1RfAfIRfBfIRfWUfIR,MfIRfVfIR4fCOCfIRfRfNNfNMfNLfO 

1 *. FALSE. , UK fIRfSSfIRfKKfIRfGfIR,P,IRfF,EIG1fUKCf JNfJXNf JUN) 
CALL  UGETIO( 3f0  f I TTY) 

507  URITE(ITTYf«)  <  'INPUT  VALUES' ) 

C  CALL  UGETIO( 1 ,03,33) 

CALL  USWFM < 4HR0EE « 4 1 ROEE ,  1  f  1 , 1  f  4  > 

CALL  USWFM ( 4HR0EC ,  4  ,  ROEC , 1 , 1 , 1 , 4 ) 

CALL  USWFM (8HMATR IX  AfBf A , IRfNNf NN , 4 ) 

CALL  USWFM <9HMATR IX  DAf9fDAf IR,NN,NNf4 » 

CALL  USWFMOHMATRIX  BfSf  B  f  IR  fNNfNMf  4  ) 

508  CALL  USWFM <9HMATR IX  DBf9fD8f IRfNNfNMf4> 

CALL  USWFM (3HMATR IX  C f B fC f IR f NK fNN f 4 ) 

CALL  USWFM (9HMATR IX  DC f9 t DCf IR ,NK fNN f 4 > 

CALL  USWFM I8HMATR IX  DfBfDf IRfNN fNQf 4 ) 

CALL  USWFM (9HMATRIX  DDf9f DDfIRf NNf NQf 4 > 

CALL  USWFM (8HMATR IX  MfBfMf IRfNLfNNf4 ) 

CALL  USWFM (9HMATRIX  DMf9fDMfIRfNLfNNf4> 

CALL  USWFM (8HMATRIX  WfBf W. IRfNQfNQ f4 > 

CALL  USWFM (8HMATRIX  VfBfVf IRfNLfNL f4 ) 

CALL  USWFM (8HMATRIX  QfBfQf IRfNKfNKf4) 

CALL  USWFM <8HMATR IX  Rf8fRfIRfNMfNMf4) 

*  URITECITTYf*) < 'OUTPUT  VALUES') 

CALL  USWFM (8HMATR IX  GfBfGf IRf NMfNN f4 ) 

CALL  USWFM (8HMATR IX  FfBfFfIRfNN.NLf4> 

CALL  USWCM< 16HEIGEN  VALUES/OPL f 16fEIG2f IR4fNNf 1 ,4 ) 

CALL  USUCM< 16HEIGEN  VALUES/NCL f 16, EIG1 f IR4 ,NN2 f 1 f4 ) 

C 

C  NOMINAL  AND  PERTURBED  CLOSED  LOOP  COSTS 

C 

10  CALL  MULT  (B.0.A12f NNfNMfNNf IRf IRf IR) 

IS  CALL  MULT<FfMfA21fNNfNLfNNfIRfIRfIR> 

20  CALL  NADDSB<Af A12fA21 fNNf A22f IRr IRf IRr IR) 

23  CALL  MULT(DB*GfA32f NNfNMfNNf IRfIRfIR) 

30  CALL  MAT  ADD  <  A  f  D A , A33  f  NN  f  NN  f I R ) 

35  CALL  HATADD(A12fA32fA34fNNfNNfIR> 

40  CALL  MULT<F fDMf A41 fNNfNLfNNf IRfIRfIR) 

45  CALL  MATADD<A21fA41fA43fNNfNNfIR> 

CALL  MATFT2(AfA12fA21fA22fNNfNNfNNfNNfAF11fNN2fNN2fIR2fIRfIR> 
CALL  MATFT2<DAfA32fA41fZER0fNNfNNfNNfNNfAF21fNN2fNN2fIR2fIRfIR> 
CALL  MATFT2<A33fA34fA43fA22fNNfNNfNNfNNfAF22fNN2fNN2fIR2fIRfIR> 
52  CALL  MATFT2<DfZER0fZER0fFfNNfNNfNQfNLfDF1fNN2fN0L 

1 fIR2fIRfIR) 

60  CALLMATFT2 ( DD  f  ZERO  f  ZERO  f  ZERO  f  NN  r NN  f  NO  f  NL  f  DF2  f NN2 

1 fNQLf IR2f IRf IR) 

WRITE (Sf*) <  <AF1I ( If J) f J*I fNN2>  f I* I fNN2) 

WRITECSf* ) ( ( AF21 (IfJ)fJ"1fNN2)»I*1 fNN2> 

65  CALL  MATADD(CfDCfC23fNKfNNf IR) 

70  CALL  MATFT4 ( C  f  ZERO  f  ZERO  f ZERO  f  DC  f  ZERO  f  C23  f  ZERO  f  ZERO , G  f  ZERO  t ZERO 
IfZEROfZEROf ZERO fGfCFfNKfNKfNMfNMfNNfNNfNNfNNfNKM2fNN4f IR, IR 
1 f IRf IRf IR4) 

90  .CALL  MATFT2<Wf  ZERO , ZERO, VfNQ f NL  fNQfNQ  f WF fNOL fNQLf IR2f IR , IR) 

C  SOLVE  LIAPUNOV'S  EQUATION! 

C  <  A)  PF1I.AF11T+AFI1.PF1I  +  DFU.U.DF11T-0 

95  CALL  MP32(IR2fNN2fDF1fIR2fN0LfWFfIR2fN0LfNN2fDF1 

1 , IR2fPR1 ) 

CALL  MEQ<AFUfDUMfNN2fNN2f  1R2,  IR2) 

100  CALL  EIGRF<DUMfNN2fIR2f1fEIG1fM0D1fIR4fWKfIER> 

105  CALL  CINV(MOD1fIMOD1fWKfIR4fNN2> 

110  CALL  LYAP2<EIG1fIR4fMODIfIMODIfIR2fPR1fIR2fPFU 

1 ,NN2fWKC ) 


C  U<)  i  h  2.:T+AF11.FF12H  K 1 1 . AF21T+&F 1 . W . HF2T-0 

115  CALL  MULRRT<PFll.AF2l.PR21.NN2.NN2,NN2.IR2.IR2.IR2) 

120  CALL  MP32  < XR2.NN2.DF1 . IR2.NQL.WF, IR2,NQL»NN2, DF2, XR2 

1.PR22) 

125  CALL  MATADD<PR21.PR22,PR2.NN2.NN2.IR2> 

CALL  MEO < AF 22 1 DUH . NN2 »NN2» IR2. IR2) 

130  CALL  EIGRF(DUM,NN2, IR2, 1 .EXG2.M0D2, IR4.UK, IER > 

135  CALL  CINV<M0D2.IM0D2.WK,IR4,NN2> 

140  CALL  LYAP5(NN2,EIG1, IRA. MODI . IM0D1 .NN2.EIG2. IR4.M0D2 

1 » IM0D2 . IR2. PR2 1 IR2  f PF 12  f  WKC  > 

CALL  USWCM< 16HEIGEN  VALUES/PCL. 14.EIG1 . !R4r NN2. 1 »4) 

C  <C)  PF22« AF22T+AF22.PF22+TPF12. AF21T+AF21 .PF12 

C  +DF2.U.DF2T  -  0 

DO  145  1“1 ,NN2 
tO  145  J«*l  ,  NN2 
145  TPF12< I , J)»PF12< J» I ) 

150  CALL  MULRRT(TPF12,AF21,PR31,NN2»NN2.NN2. IR2» 1.12# IR2) 

155  CALL  MULT<AF21.PF12.PR32.NN2»NN2.NN2,IR2.IR2.IR2> 

140  CALL  MP32<IR2,NN2,DF2.IR2,N0L,UF,IR2.NaL.NN2»DF2.IR2 

1 t PR33) 

145  CALL  MATAtD<PR31»PR32,DUM,NN2.NN2, IR2> 

170  CALL  MATADD<tUM,PR33.PR3.NN2.NN2,IR2> 

185  CALL  LYAP2<EIG2»IR4.M0D2.IM0D2,IR2,PR3,IR2.PF22,NN2,UKC> 

C  FORM  MAT  PFJ 

195  CALL  HATFT2<PF11,PF12,TPF12.PF22.NN2»NN2.NN2.NN2.PF»NN4 

lfNN4fIR4.IR2.IR2) 

C  FORM  CFT.OFl.CFf  CFT.0F2.CF.  CFT.GF.CF 

194  CALL  MATFT4 <0,0. ZERO . ZERO . Q , 0. ZERO . ZERO . ZERO , ZERO . ZERO . ZERO . ZERO . ZERO 
1, ZERO. ZERO. OF .NK.NK»NH»NM.NK»NK.NH,NM,NKM2»NKM2 , IR » IR, IR , IR , IR4 > 

200  CALL  MP31 ( IR4 . NN4 . CF » IR4 .NKM2.QF . IR4,NKM2. NN4 ,CF . IR4 . COC ) 

201  JXP-TRAC2  < IR4.PF , IR4, COC , NN4 . NN4  > 

202  CALL  MATFT4 <  ZERO . ZERO . ZERO , ZERO . ZERO . ZERO . ZERO » ZERO . ZERO . ZERO 
1. R.R. ZERO. ZERO, R,R. OF. NK.MK.NM.NM.NK.NK.NM.NM.NKM2.NKM2. IR, IR 
1 , IR , IR , IR4 ) 

205  CALL  MP31<IR4,NN4,CFf IR4.NKM2.QF. IR4.NKM2.NN4.CF. IR4.C0C) 

204  JUP-TRAC2  < IR4.PF, IR4 ,CQC  f NN4 .NN4 ) 

207  CALL  MATFT 4 (0,0, ZERO , ZERO ,0,0. ZERO , ZERO , ZERO , ZERO , R.R, ZERO , ZERO , 

1R.R, OF, NK.NK »NM,NM, NK.NK .NM.NM.NKM2 »NKM2. IR. IR, IR, IR. IR4) 

210  CALL  MP31 ( IR4, NN4, CF , IR4 .NKM2.QF. IR4 » NKM2 , NN4 , CF , IR4, COC) 

225  JP*  TRAC2  < I R4 , PF  » IR4 , COC , NN4 , NN4  > 

CALL  MATFT4 ( 0 , ZERO , ZERO , ZERO , ZERO , ZERO , ZERO , ZERO , ZERO , ZERO 
1,R, ZERO, ZERO, ZERO. ZERO. ZERO. OF, NK.NK. NM,NM,NK,NK.NM,NH,NKM2 
1 ,NKM2, IR.IR.IR, IR, IR4 ) 

1001  CALL  MP31  <  IR4,NN4»CF  .IR4.NKM2.QF  .1R4.NKM2.NN4  ,CF  ,  IR4.CQC) 

1002  JNN-TRAC2 ( IR4 ,PF  » IR4 , CQC , NN4 , NN4  > 

1003  CALL  MATFT4(0, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO 
1.  ZERO,  ZERO,  ZERO.  ZERO,  ZERO.  ZERO.  OF,  NK, NK.NM, NM.NK.NK.NM.NM 
1 , NKM2 . NKM2 ,IR,IR,IR,IR,IR4> 

1004  CALL  MP31(IR4,NN4,CF , IR4.NKM2.0F , IR4.NKM2.NN4 ,CF,IR4,CQC) 

1005  JXMN-TRAC2 (IR4,PF,IR4, COC , NN4 , NN4 ) 

1004  CALL  MATFT4< ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO 
1.  ZERO.  R.  ZERO,  ZERO,  ZERO.  ZERO,  ZERO.  QF»NK,NK,NM,NM,NK.  NK.NM, 
1NH,NKM2,NKM2,IR, IR, XR, IR, IR4) 

1007  CALL  MP31 ( IR4.NN4.CF , IR4.NKM2.QF, XR4.NKM2, NN4 ,CF , IR4.CQC) 
JUNN-TRAC2(IR4,PF,IR4.CQC,NN4,NN4) 

BETAPR-ABS < JXP- JXNM ) / JXNN 
CALL  XDENT(Q,NK,IR> 

1008  CALL  MATFT4  (  0 ,  ZERO  ,  ZERO ,  ZERO .  ZERO .  ZERO  ,  ZERO  ,  ZERO .  ZERO ,  ZERO ,  ZERO 
1 ,  ZERO  ,  ZERO .  ZERO ,  ZERO  ,  ZERO ,  OF ,  NK ,  NK  ,  MM ,  NM  ,  NK  ,  NK ,  NM  ,  NM  ,  NKM2  ,  NKM2 , 

HR, IR.IR.IR, IR4) 

CALL  HP31 < IRA , NN4 , CF , IRA , NKM2 , OF , IRA  ,NKH2 , NN4 , CF , XR4 , CQC ) 

1009  YTY-TRAC2 ( IR4 , PF , IRA , CQC , NN4 , NN4  > 

1010  CALL  MATFT4( ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO 

1, ZERO, Q, ZERO, ZERO, ZERO, ZERO, ZERO, OF, NK, NK.NM, NM.NK.NK.NM.NM 
1 ,NKM2,NKM2, IR.IR.IR, IR, IRA) 

1011  CALL  HP31 (IR4,NN4,CF , IR4.NKM2.0F » IR4.NKM2, NN4 »CF , IR4.CQC ) 

UTU-TRAC2 < IRA , PF , IRA . CQC » NN4 . NN4 ) 

MRITEdTTY.tX'OPEN  LOOP  COST') 

CALL  USWFM< 14H0PEN  LOOP  COST, 14, JOP, 1 , 1 , 1 , A) 

WRITE! ITTY»*>  < 'NOM. CLOSED  LOOP  COSTS') 

CALL  USUFM( 15HRE0ULATI0N  COST , 15, JXNN, 1 , 1 , 1 , A) 

CALL  USUFMf 12HC0NTR0L  COST, 12. JUNN.l , 1 , 1 ,4) 

CALL  USUF?1<10HTOTAL  COST,  10, JNN,  1 , 1 ,1 ,4 ) 

WRITE(ITTY.*>< 'PERT. CLOSED  LOOP  COSTS') 

CALL  USWFMU3HREBULATI0N  COST,  13, JXP.l ,  1 , 1 . A) 

CALL  USMFN( 12HC0NTR0L  COST, 12, JUP, 1 ,1 , 1 , A) 

CALL  USWFM< 10HT0TAL  COST, 10, JP, 1 , 1 , 1 ,4) 

CALL  USWFM(3HY  T  Y,5, YTY, 1 , 1 ,1 ,4) 

CALL  USWFM(5HU  T  U.3.UTU, 1 , 1 . 1 ,4) 

TTY1«SQRT<YTY) 

UTU1 "SORT ( UTU ) 

CALL  USWFM  ( 8HYTYM1/2 . 8  ,  YTY1 , 1 , 1 , 1 , 4 ) 

CALL  USUFM  <  8HUTUM 1  /2 , 8  .  UTU1 . 1 . 1 . 1 »  4  > 

MR I TE ( I TT Y . * > < ' PERF . ROB . I NDEX ' > 

CALL  USWFM(7HBETA  PR,7,BETAPR,1 » 1 ,1 ,4) 

FORMAT (F18. 5) 

FORMAT (3F18. 5) 

STOP 


12 

13 


36 


SUBROUTINE  HADDSB  t A»Br C .N, ABC. IA. IB. ICr IABC) 

C  THIS  DOES  ABC-A+B-C 

REAL  AdA»N)»BdB»N)rC( ICtN)  •  ABCdABC>N> 

DO  41  I  ■  liN 
DO  41  J  Mi  N 

41  ABC<If J)-A<X'J)+B<Xr J)-C(ItJ) 

RETURN 

END 


SUBROUTINE  MATFIT  < A1 1 t A12* A13 r A21 * A22 r A23. A31 * A32 r A33 
1  *NRl.NR2»NR3»NClfNC2*NC3»NR.NC»Af IR1»IR2» IR3.  IR> 

REAL  All(XRlrNCl>'A12<lRl'NC2>.A13(XRl'NC3>« 

1  A21 <  XR2>NC1 ) p  A22dR2r  NC2) > 

1  A23<  XR2rNC3> pA3ldR3»NCl ) »A32< IR3# NC2) » A33( XR3.NC3) • A< IR .NC) 

DO  40  1-1 »NR 
BO  40  J-1»NC 

40  A<I»J)>0 

C  THIS  SUBR  COMPILES  MATRICES  INTO  A  ONE  MATRIX 
DO  .41  I  -lrNRl 
DO  41  J  *1 t NCI 
A(IfJ)  -  All  <I»J) 

41  CONTINUE 

DO  42  I  -l'NRl 
DO  42  J  -1.NC2 
A(I»NC1+J)-A12<I»J) 

42  CONTINUE 

DO  43  I  *1 pNRI 
DO  43  J  -1.NC3 

43  Ad»NCl+NC2+J>«A13dpJ> 

DO  44  I  -1pNR2 

DO  44  J  «1»NC1 

44  A(NR14I»J>  -  A21 < I » J) 

DO  45  I  »1»NR2 

DO  45  J  -1pNC2 

45  A<NR1+I*NC1+J>-  A22(I»J> 

DO  44  I  -1.NR2 

DO  44  J  -1.NC3 

44  A<NRl+IfNCl+NC2+J)-A23d»J> 

DO  47  I  -1.NR3 
DO  47  J  -If NCI 

47  A<NRl+NR2+X'J)-A31df J) 

DO  48  I  -  I.NR3 

DO  48  J  -Ip  NC2 

48  A<NR1+NR2+I»NC1+J>-A32dr  J> 

DO  49  I  ■  1 »NR3 

DO  49  J  NC3 

49  A(NRl+NR2+I»NCl+NC2+J>-A33d*J> 

RETURN 

END 


8UDR0UTXNE  ABS* (AApNRpNCpIA) 

ctmmnmtmmmmmmmmmmttmmmm 

C  THIS  8IVES  ABSOLUTE  VALUES  OF  ALL  ELEMENTS  OF 
C  AA  MATRIX  AND  OUTPUT  IS  STORED  BACK  INTO  'AA' 

REAL  AA<XA?NC) 

DO  71  1*1 * NR 
DO  71  J>1»NC 

71  AAIXr  J)«ABS(AA(  X#  J) ) 

RETURN 

END 

SUBROUTINE  DXAMAT<AD?NDfEXGrIADpIEXG) 

•  ■  ,**«*« 

C  THIS  SUBROUTINE  HAKES  'AD'  A  DIAGONAL  MATRIX 
C  DIAGONAL  ELEMENTS  AS  REAL  PART  EIGEN  VALUES  OF 
C  ANOTHER  MATRIX  OF  SIMILAR  ORDER 
REAL  ADdADtND) 

COMPLEX  EXGdEXG) 

DO  SI  I«1 *ND 
DO  81  J«lrND 

81  ADd>J)>0 

DO  82  I-1pND 

82  ADd.I)  ■REALIEIGd) ) 

RETURN 

END 


37 


SUBROUTINE  MATFT4 


C 

c 

c 

C  *t* 

c  ««> 


c  **« 

C  Mt 


THIS  SUBROUTINE  COMPILES  4*4  MATRICES  INTO  A  SINGLE  MATRIX. 


c 

•*« 

All 

A12 

A13 

A14 

c 

»*« 

A21 

A22 

A23 

A24 

c 

III 

A31 

A32 

A33 

A34 

c 

Ml 

A41 

A42 

A43 

A44 

C  Ml 
C  *** 
C  Mt 
C  *** 
C  *«* 
C  *•* 

c  *«« 
c  «** 
c  *** 
c  »*** 
c 


A14-NR17NC4 

A24-NR24NC4 


All -NR I INC 1 »  A12-NR1INC2*  A13-NR1INC3. 

A21-NR2INC1 »  A22-NR2INC2.  A23-NR2*NC3» 

A31-NR3INC1.  A32-NR3INC2 t  A33-NR37NC3,  A34-NR3INC4 
A41-NR4INC1.  A42-NR4INC2 »  A43-NR44NC3,  A44-NR4INC4 
A-NRINC 


MR1 #MR2»MR3f MR4 1 MCI r hC2  »MC3>  MC4 ,HR > MC-DIMENSIONS 


SUBROUTINE  MATFT4< All t A12f A13>A14>A21 • A22»A23'A24>A31 t A32 

1 1  A33>A34>A41  r A42»A43r A44«ArNRl  r NR2 r NR3 >NR4 r NCI  >NC2rNC3rNC4 
1 »NR  t NC  >MR1 »MR2  >  MR3 *  MR4  » MR) 


REAL  All <MR1 t MR1 ) r A12 (MR1 »MR2 >»A13(MR1 »MR3> >A14(MR1>MR4) 
REAL  A21 < MR2 »  MR1 > . A22  <  MR2 1 MR2 ) » A23  < MR2 » MR3 ) . A24 ( MR2 » MR4  » 
REAL  A31 < MR3 * MR1 ) . A32 <MR3» MR2 > t A33 < MR3 . MR3 ) , A34 ( MR3. MR4 ) 
REAL  A41 <  MR4  *  MR1 ) t A42 < MR4 • MR2 >  » A43IMR4  »MR3) »A44(MR4,MR4) 
REAL  A (MR t MR) 

C 

DO  SO  I-l.NR 
DO  SO  J»1>NC 
SO  A< I t J)«0 

DO  100  I>1.NR1 
DO  100  J-l.NCl 
100  A(I>J)»A11(I»J) 

C 

DO  110  X-lrNRl 
DO  110  J-1.NC2 
110  A<X,NC1+J)-A12U'J) 

c 

DO  120  I*1»NR1 
DO  120  J-1.NC3 

120  A<I>NCl+NC2+J>-A13<Ir J) 

C 

DO  130  X»1.NR1 
DO  130  J-1.NC4 

130  A(I>NCl+NC2+NC3+J)-A14(Ir,J> 

c 

DO  140  X-lrNR2 
DO  140  J*lfNCl 
140  A(NR14I»J)«A21(Xf J) 

C 

DO  ISO  I-lrNR2 
DO  ISO  J-1.NC2 

150  A(NRl+I,NCl+J)*A22(IfJ)  ' 

C 

DO  140  X»1 ?NR2 
DO  140  J-1.NC3 

if  M#»4  •  J  .  |\ 

DO  170  I>1»NR2 
DO  170  J»1»NC4 

170  A<NR1+I,NC1+NC24NC3+J)-A24<I.J> 

c 

DO  180  X-1.NR3 
DO  180  J-1,NC1 

180  A<NR14NR2+IrJ)-A31<X,J> 

C 

DO  170  I-1.NR3 
DO  170  J"1 t NC2 

170  A<NR1+NR24X»NC14J)-A32<I.J> 

C 

DO  200  X-lrNR3 
DO  200  J»ltNC3 

200  A  <  NR 1 +NR2+ 1 , NC 1 4NC2+ J ) • A33  <  X  »  J  > 

c 

DO  210  I-1.NR3 
DO  210  J-1.NC4 

210  A(NR1+NR24I r NC1+NC24NC3+J)*A34(X  r J> 
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00  220  J-1.NC1 

220  ACNR1+NR2+NR3+I • J)-A41  ( I  >  J) 

C 

00  230  X-1.NR4 
00  230  J-1.NC2 

230  A<NR1+NR2+NR3+I .NC1  + J) «A42< I . J> 

C 

00  240  1*1 .NR4 
DO  240  J*1 »NC3 

240  A<NR1+NR24NR3+I  *NC1+NC2+J)*A43<  1  *  J> 

C 

DO  2S0  I-1.NR4 
DO  230  J-1.NC4 

230  A(NR1+HR2+NR34’I  .NC1+NC2+NC3+ J)-A44<  I »  J) 

C 

RETURN 

END 

SUBROUTINE  NATFT2  ( All , A12 . A21 r A22, NR1 .NR2 r 
1NC1 .NC2. Ar NR.NCr IR. IR1 »  XR2) 

C  THIS  FORMS  A  MATRIX  FROM  All .A12rA21 . A22 

C 

C  All  A12 

C  A21  A22  TO  FORM  A 

C  NR-NR1+NR2I  NONC1+NC2I 
C 

C  NR-NR1+NR2 

C  NC-NC1+NC2 

C 

REAL  A(IR.IR).A11<XR1.XR1).A12(XR1.XR2) 
REAL  A2KXR2.1R1)  . A22( XR2. XR2) 

DO  10  1*1. NR 
DO  10  J*1.NC 
10  A<I.J)-0 

DO  1  I-1.NR1 

DO  1  J— 1.NC1 

1  A(X.J)  -All( I . J) 

DO  2  1*1. NR1 

DO  2  J-1.NC2 

2  A(X.NC1+J)*A12<I.J) 

DO  3  Z-1.NR2 

DO  3  J-l.NCl 

3  A(NR1+I . J)-A21 ( I. J) 

DO  4  X-1.NR2 

DO  4  J-1.NC2 

4  A<NR1+I .NC1+J)*A22< I > J) 

RETURN 

END 


SUBROUTINE  STABLER A. B.C.HK.NU.D.AE.N.U.Z 
1  . 1A.IB.IC.ID.IAE.IZ) 

ct*tt***************s***t*t**t*******«****t*t****t**t**t*tt**t* 

C  THIS  SUB  FINDS  EIGEN  VALUES  OF 
C  C  (E)St(INU(AD) >S  3SYM 

REAL  A(IA.N).B(IB.N).C(IC.N).UK(NU).D(ID.N) »AE ( IAE.N) 

INTEGER  N.IJOB.IER 

COMPLEX  W<IZ).Z<XZ.N)»ZN 

X JOB-O 

DO  SS  X-l.N 

DO  S3  J-l.N 

SS  AE<I.J)-(A(X.J)+A< J.X) )/2 

C  FORMS  XNV  DIAGONAL  MATRIX  ONLY 
DO  44  1*1, N 
DO  66  J-l.N 
IF(I.EQ.J)  GO  TO  44 
B(I. J)-B(I. J) 

GO  TO  44 
44  CONTINUE 

B<X.X)-1«0/B(X.X) 

44  CONTINUE 

CC  HRITE<S.*)('AD  MATRIX') 

C  URITE(5.*)((B<I.J). J-l.N). I-l.N) 

CALL  VMULFF  (AE.B.N.N.N. IAE. IB.C. IC» IER) 

C  WRITE (5. *> ( 'AE  MATR ' ) 

C  WRITE <3. t) <(AE(I.J) » J-l.N). I-l.N) 

C  WRITE<5.*K'C  MAT') 

C  WRITE <5. * ) < (C(I.J). J-l.N). X-l.N) 

DO  9  I-l.N 
DO  9  J-l.N 

9  D<I»J>-  <<C<X.J)+C(J.I))/2) 

CALL  CXORF  (D.N.ID.IJOB.W.Z.XZ.WK.IER) 

END 
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SUBROUTINE  MARS<EIO» IR»N»SMALL > 

C  THIS  SUBROUTINE  FINOS  MIN  VALUE  REAL  PART  OF  COMPLEX 
C  VALUES 

C  REAL  REL(IR) 

COMPLEX  EIQ< IR) 

SMALL-REAL < EIG ( 1 )  ) 

6MALL-ABS< SMALL) 

DO  1  I-l.N 
REL-REAL<EIG<I)) 

REL-ABS(REL) 

IF (REL.GE. SMALL)  GO  TO  1 
SMALL-REL 
X  CONTINUE 
RETURN 
END 


SUBROUTINE  MAR8<EIGf IRtN.BIG) 
COMPLEX  EIG(XR) 

BIGI-REAL<EIG< 1 ) ) 
BIG-ABS(BIGI) 

DO  1  I-1fN 
REL-REAL<EIG<I>) 

REL»ABS(REL) 

IF (REL.LE.BIG)  GO  TO  1 
BIG-REL 

X  CONTINUE 
RETURN 


SUBROUTINE  8TABR<ACLfECLfALPHAfIDfAEfEXGfZfNf 
IXRfXZfIEXGfUKfIUK) 

REAL  ACL(XRfN)fECL(XRfN) fID(XRfN) fAE(IRfN) rUK(IUK) 
COMPLEX  EIG(IEIG) rZ(IZrN) rZH 
I JOB-2 

12  CALL  XDENTUDfNfXR) 

13  CALL  SCAHUL(ALPHAfXRfIDfXRfXDfNfN) 

14  CALL  MADDIACLfIDfAEfNfNfIRfIRfIRfI) 

13  CALL  ABS(AEfNfNfIR) 

16  CALL  MADDIAEfECLfAEfNfNfXRfXRf IRfI  ) 

DO  1  X-1*N 

DO  1  J-l.N 

1  ID( I» J>-< AE( Ir J)+AE( J» I ) )/2 

ALP HAI-1 /ALPHA 

17  CALL  8CAMUL( ALPHA! f  XRfXDfIRf IDfNfN) 

15  CALL  EIORF(XD.NFXRFXJOBFEXGFZFXZrUKFXER) 

RETURN 

END 

SUBROUTINE  MAR<EXG«XRfNf SMALL) 

C  THIS  SUBROUTINE  FINDS  NXN  VALUE  REAL  PART  OF  COMPLEX 
C  VALUES 

C  REAL  RELIIR) 

COMPLEX  CXO(XR) 

SMALL-REAL ( EIG < 1 ) ) 

DO  1  I-l.N 
REL-REAL(EXGd)  ) 

XriREL.OE. SMALL)  GO  TO  1 
SMALL-REL 
CONTINUE 
END 


1 


PROGRAM  TEST 

REAL  ACL<0t8) rECL(8>8) tCC<S»G)»AE(8»8) f UK <1000) 
COMPLEX  EICN(B) >EE (8 >8) 

18-8 

N-4 

WRITEtSfBX'GIVE  SIZE  OF  MATRIX') 

READ<3f«)N 

NRIT£<3f*X 'INPUT  VALUES  OF  ACL') 
R£AD<36f*X<ACL(If JXJ-l.N) ,1-1. N> 

WRITE<5f*X 'INPUT  VALUES  OF  ECL'> 

REA0(36>  t ) ( (ECL< IfJ)rJ»lfN)rI-l>N) 

DO  1  I  -1  f 23 

WRITE<5f*X'GIVE  VALUE  OF  ALPHA') 

READ <  3  f  *  > ALPHA 

CALL  STABR ( ACL f ECL » ALPHA .CC.AE>EIGN>EE>N> 

IXBf IBf IBfUKf 1800) 

CALL  UGETI0(3f Ot 5) 

•  CALL  USUFM(3HACLt3fACLf I8fNfN>4) 

CALL  USUFM(3HECLt 3>ECLf IB>Nf Nf 4 ) 

CALL  USUCM  <£HEIGENSf6f EIGNr IBf Nf 1 f 4) 

CALL  MARB<£IGNf IRf Nf BIG) 

URITE(3f*X '/MAX  EIGEN/  -  ' f BIG) 

IFtBIGfLT# l)GO  TO  10 

IF(BIG.GE.l)URITE<Sf*X'SYS  IS  NOT  STABLE') 
WRITE(3f <X 'TYPE  1  TO  CONTINUE') 

HRlTEISf B) ( 'TYPE  2  TO  STOP') 

READ(3f*)NT 
IF(NT.EO.l)  GO  TO  1 
GO  TO  100 
CONTINUE 
CONTINUE 

URITE(Sf*X'SYS  IS  STABLE') 

CALL  EIGRFCACLrNf IBt 2f EIGNf EEf IBt UKf IER) 

CALL  USWCM<6HEXGACLf AfEIGNf IBfNrl f4) 

CALL  MARS(EIGNf IBfNf SMACL) 

■  WRITE<3f«X 'SMACL  -'f SMACL) 

CALL  HATADDIECLf ACL.ACLfNf Nf IB) 

CALL  EXGRF(ACLfNfXBf2fEIGNfEEf IBfUKflER) 

CALL  USUCMCAHACLECLf AfEIGN. IBfNf If 4) 

CALL  MARSIEIGNf IBfNf SMPER) 

WttTE(3«tU'6MPER  «' .SUPER) 

BSTAB- <  ABS ( SMACL-SMPER ) ) /SMACL 
URITE(3f 8) ( 'BSTAB*  'f BSTAB) 

CONTINUE 

STOP 


)lication  Example: 


Consider  the  Scalar  System 

00 

2  2 

x  =  -  x  +  u,  x  =1,  J  =  (x  +pu)dt 
’on  c  ' 

0 

i.e.  a  =  -1. ,  b  =  1 

Let  |Aa|  =  0.5;  (Abj  =  1.207.  The  closed  loop  system  is  given  by 

x  =  (a  +  bg)  x  =  aCL  x 

and  the  perturbation  is  given  by  Aa^. 

We  now  consider  three  cases  of  the  perturbation  Aa^- 

Case  I:  •  AaQL  =  |Aa|  ♦  |Ab|g 

Case  II:  Aa£L  =  | Aa |  +  [ Ab |  |g| 

Case  III:  AaCL  *  -  (|Aa|  +  |Ab|  |g  [) 


The  robustness  indices  8C  D  and  SD  D  are  calculated  for  different 

values  of  p  and  tabulated  in  Table  I.  Note  that  all  the  controllers 
c 

corresponding  to  the  range  of  pc  considered  satisfy  the  stability  ro¬ 
bustness  condition. 

Note  that  the  testing  of  stability  condition  is  done  with  the  per¬ 


turbation  A^  corresponding  to  Case  II  to  satisfy  theorem  B.  Then  by 
virtue  of  theorem  B,  if  the  perturbed  system  with  Case  II  perturbation  is 
stable  then  the  perturbed  systems  of  Cases  I  §  III  will  also  be  stable. 


_  *5  T  H  Case  I  Case  II  Case  III 

c  (X  x)  U  U  gS.R.  SP.R.  ^S.R.  ^P.R.  0S.R.  ^P.R. 


■ 


.594  10.246 


0.39 

0.97 

0.31 

0.89 

0.18 

0.798 

0.08 

0.744 

0 

0.707 

0.07 

0.680 

0.132 

0.658 

1/4  0.473  10.584  0.443  0.31  0.89  9.15  0.89  0.47 


1/2  0.537  0.393  0.22  0.18  0.798  3.97  0.798  0.44 


3/4  |o .572  0.302  0.08  I  0.08  0.744  2.91  0.744  0.43 


5/4  0.610  0.209  0.065  0.07  0.680  2.125  0.680  0.405 


3/2  jo. 6225  0.181  0.116  0.132  0 .65S|  1.93  0.658  0.397 


2  10.639  0.143  0.187 


BE 


1.70  10.63  0.385 


From  the  above  table,  it  can  be  observed  that  Case  II  gives  the  worst  case 
8C  _  and  worst  case  $D  „  .  Thus,  we  use  these  indices  for  comparison  of 
different  controllers  from  stability  and  performance  point  of  views.  The 
information  corresponding  to  Case  II  is  presented  graphically  in  figures 
A,B,C  and  0.  A  reasonable  choice  for  robust  controller  could  be  the  one 
corresponding  to  pc  =  1/2,  3/4  or  1. 

Since  this  is  a  scalar  example  it  is  easy  to  see  which  case  of  perturbation 
represents  the  worst  case  situation.  This  may  not  be  easily  inferred  for 
matrix  cases.  To  determine  which  case  of  perturbation  (among  the  possible 
three  cases  considered)  represents  the  worst  case  situation  for  general 
matrix  case  is  suggested  as  a  future  research  topic. 

-  43  - 


III.  Conclusions  and  Recommendations  for  Future  Research: 


The  main  theme  of  the  Mini  grant  research  has  been  to  investigate  further 
into  the  development  of  a  stability  robustness  condition  in  time  domain  and 
the  extension  of  these  results  in  the  computer  implementation  of  a  robust 
control  design  algorithm  that  incorporates  both  stability  robustness  and 
performance  robustness  into  the  control  design  procedure.  Towards  this 
direction,  first  a  new  stability  robustness  condition  is  developed  in  time 
domain  (in  terms  of  eigenvalues)  and  it  is  shown  that  the  proposed  time 
domain  condition  is  less  conservative  than  the  corresponding  frequency 
domain  condition  as  well  as  another  recently  developed  time  domain  con¬ 
dition.  Also,  further  observations  are  made  in  the  comparison  of  proposed 
time  domain  development  to  the  frequency  domain  development.  Then  new 
measures  of  'stability  robustness'  and  'performance  robustness'  are  developed 
and  incorporated  into  the  robust  control  design  algorithm  proposed  in  the 
summer  research.  Finally,  computer  software  is  developed  to  implement  the 
proposed  control  design  algorithm  and  examples  are  presented  which  involve 
the  use  of  the  software. 

The  experience  with  Large  Space  Structure  examples  carried  out  in¬ 
dicates  that  for  these  models  the  stability  condition  in  its  present  form 
is  still  conservative  and  that  more  research  is  needed  to  specialize  the 
analytical  development  to  LSS  models. 


As  it  normally  occurs,  another  result  of  this  study  is  that  many  in¬ 
teresting  research  topics  surfaced  for  further  investigation.  In  the  next 


X  '.N  .I.  .*.  .S  .V 
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section,  we  project  those  areas  which  merit  serious  research  effort  in  the 
form  of  shponsorship  from  AFOSR. 

Avenues  for  Further  Research: 

1)  The  foremost  area  of  research  would  be  to  look  further  into  the  ' 
stability  robustness  condition  from  two  viewpoints. 

a)  To  improve  the  ’optimism*  of  the  proposed  condition  particularly 
with  reference  to  LSS  models. 

b)  To  investigate  the  possibility  of  developing  a  new  stability 
robustness  condition  which  is  both  necessary  and  sufficient. 
(Recall  that  the  present  condition  is  a  sufficient  condition) . 

2)  Another  area  of  immediate  concern  is  to  arrive  at  an  algorithm 

(a  technique)  that  would  give  the  worst  case  (Jg  R  and  worst  case 
3d  for  given  perturbations,  for  comparison  purposes. 

3)  An  area  of  extreme  interest  would  be  to  develop  an  algorithm  for 
’Maximum  Allowable  Perturbations'  that  would  destabilize  a  given 
stable  system.  In  a  way  this  is  an  ’inverse’  problem.  This 
problem  is  apparently  related  to  the  task  number  one  indicated 
above. 

4)  An  important  area  of  research  is  to  investigate  any  computation 
reduction  schemes  for  the  proposed  algorithm. 

5)  It  is  also  of  interest  to  probe  further  into  the  relationship 
between  frequency  domain  treatment  and  the  proposed  time  domain 


treatment . 
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Proof  of  Theorem  1 : 


Let  p{(Es(Fs)*1]s)<l 

-  |A(Es(Fs)  _1)J  <  1 

max 

-1 

-  | A.  (E  (F  )  )  |  <  1  Vi 

-  1  +  >  0  Vi 

-  A.{l  +  (Es(Fs)_1)s}  >  0  Vi 

-1 

-  A.{[I  ♦  Es(Fg)  Js}  >  0  Vi 

[I  +  E  (F  )  *]  is  positive  definite  ** 

s  s 

-*■  [I  +  E  (F  )  X]  I-F  ]  has  positive,  real  eigenvalues  because  1)  if 
s  s  s 

A  and  B  are  positive  definite,  AB  has  positive 
real  eigenvalues.  (Ref  [24,25]) 

and 

2)  If  A  is  negative  definite,  -A  is  positive 
definite  and  hence  -F  is  positive  definite 
(Ref  [26  3). 

■*  -  (Fg  +  Es)  has  positive,  real  eigenvalues 

[because  [I  +  E  (F  )_1J  [-F  ]  =  -  (F  +  E  )] 

5  a  5  S  S 

-  (Fg+  Es)  is  positive  definite  (because  -  (Fg+  Eg)  is  symmetric  too) 

-*■  (Fs  ♦  Es)  is  negative  definite 

**■  (F+E)  *s  negative  definite 

-*>  (F+e)  is  dtf'nit* 

■*  (F+E)  has  negative  real  part  eigenvalues 


(F+E)  is  stable 


